以頻率誤差控制為目標的自由振動問題自適應有限元分析

孫浩涵，袁 駟

(1. 中物院高性能數(shù)值模擬軟件中心，北京 100088；2. 北京應用物理與計算數(shù)學研究所，北京 100088； 3. 清華大學 土木工程系，北京 100084)

自由振動反映結構動力特性，是沖擊等復雜力學分析的基礎。在數(shù)學模型上，結構的自由振動可歸結為偏微分方程特征值問題，通常以有限元法等數(shù)值方法進行求解，在各類實際工程分析中得到了廣泛應用[1-2]。有限元法通用靈活，但作為一種近似數(shù)值方法會引入離散誤差，網(wǎng)格劃分將決定求解的效率和精度。自適應有限元法在傳統(tǒng)有限元法(finite element method,FEM)的基礎上，引入誤差估計算法和網(wǎng)格細化技術，自適應地反復迭代和調整網(wǎng)格，最終提供優(yōu)化的有限元網(wǎng)格和滿足精度要求的解答。這一概念最早于20世紀70年代末由Babu?ka等[3-4]提出，首先被應用于線性問題，而后被推廣至各類非線性問題。

特征值問題作為一類特殊的非線性問題，通常需要迭代求解，合理的網(wǎng)格劃分將顯著減少計算量。誤差估計是特征值問題自適應分析的核心環(huán)節(jié)。在理論方面，Strang等[5]、Babu?ka等[6]、Knyazev等[7]先后給出了特征值問題的先驗誤差估計；在應用方面，更有價值的是后驗誤差估計，相關研究包括利用殘差基于對偶性構造誤差估計[8]、基于余能原理分析特征值上下界[9]、基于p型離散和Rayleigh商的Taylor展開形式構造誤差估計[10]以及基于對應線性邊值問題的誤差估計[11-12]等?；谡裥统諗科雌謴徒?，王永亮[13]提出了中厚圓柱殼自適應有限元分析策略，成功將自適應策略運用至自由振動問題。

作為一類超收斂計算方法，相較于經典SPR法[14-15]，單元能量投影法(element energy projection，EEP)不依賴超收斂點[16]，且恢復精度更高，可實現(xiàn)有限元解的逐點誤差估計。針對結構自由振動，Sun等[17]已建立了一套基于EEP法的自適應分析策略，對頻率和模態(tài)同時進行誤差控制，可連續(xù)求解若干階結果，給出滿足精度要求的頻率及逐點滿足誤差限要求的模態(tài)[18-19]。

在實際工程應用中，也存在另一類需求，即只要求頻率滿足精度，而并不關心模態(tài)誤差大小。本文提出了頻率的后驗誤差估計算法，繼而建立了整體頻率誤差和局部模態(tài)誤差的轉換關系，最終實現(xiàn)了以頻率誤差控制為目標的自適應有限元分析策略。本方法的有效性在二階Sturm-Liouville問題(簡稱SL問題)及彈性薄膜自由振動問題上得到了應用驗證。

1 背景知識

1.1 研究對象

為方便一般性的討論，本文沿用數(shù)學表達，自由振動問題中的頻率和模態(tài)分別對應于特征值問題的特征值(平方根)和特征函數(shù)，不再加以區(qū)分。問題的微分方程形式為

Lu=λru

(1)

以及相應給定的邊界條件。式中：L為定義的微分算子；λ為特征值；u為特征函數(shù)；r為質量項。

式(1)對應的能量泛函為

(2)

式中：a(·,·)為能量內積；(·,·)為常規(guī)的雙線性型。由式(2)的一階變分，得弱形式方程為

δπ(u)=a(u,δu)-λ·(u,δu)=0

(3)

1.2 單元能量投影法

作為一種有效的有限元超收斂計算算法，EEP法是本文特征值問題誤差估計的核心。其基本思想是通過假設有限元法的能量投影定理

a(u-uh,vh)=0, ?vh∈Sh

(4)

在單個單元上近似成立，推導出具有超收斂性質的EEP解u*。然而，這一性質僅對一維問題成立，對于二維乃至三維問題，應采用逐維離散、逐維恢復的方法。如圖1所示，二維問題有限元離散可視作兩次一維有限元離散過程的疊加，進而其超收斂解可通過兩次應用一維EEP技術接續(xù)計算獲得。該方法的有效性已在一系列二維及三維問題的應用中得到了確認[20-22]。

圖1 二維有限元網(wǎng)格逐維離散Fig.1 D-by-D discretization of 2D problems

2 總體策略

2.1 自適應分析目標

本文研究對象為特征值，自適應分析最終目標是：在精確解λk(k=1,2,…,n)未知的情況下，獲得充分好的有限元網(wǎng)格πk，使得有限元解滿足

(5)

式中，Tol為事先給定的誤差限。由于精確解未知，實際采用以下誤差控制準則

(6)

2.2 整體和局部誤差

在自適應有限元分析中，既需要整體誤差控制標準，以確定何時停機；又需要局部誤差控制標準，以驅動網(wǎng)格迭代。在本文中，由于整體誤差以特征值進行衡量，而局部誤差仍需以特征函數(shù)進行評估，因此需要建立關系，將二者有機地聯(lián)系起來。

采用h型網(wǎng)格加密方式的自適應求解，收斂過程往往經歷“均勻加密”到“局部加密”兩個階段，自適應收斂率約為最佳收斂率：對于一維問題，特征值收斂率為2m，特征函數(shù)收斂率為m+1；對于二維問題，特征值收斂率為m，特征函數(shù)收斂率為(m+1)/2。

(7)

式中，n(n=1,2)為問題維度。若假設新網(wǎng)格下特征值誤差剛好滿足誤差限，則由式(7)關系可得特征函數(shù)最大誤差限為

(8)

誤差大于該值的單元均應細分。在實際算法中，為保證魯棒性，誤差限總取為歷史路徑的最小值，即

(9)

2.3 整體算法

本文自適應有限元分析總體策略可總結為：

步驟1有限元解uh。在當前網(wǎng)格下(初始網(wǎng)格用戶給定)，求解有限元解。

步驟2EEP解u*。依照Gauss積分點分布確定采樣點位置，計算超收斂位移及其導數(shù)。

步驟3誤差檢驗。計算超收斂特征值，進行特征值后驗誤差估計，若滿足式(6)則停機；否則，依據(jù)式(9)確定局部特征函數(shù)誤差限，檢查所有單元，對未滿足誤差要求的單元進行標記。

步驟4單元細分。采用基于四叉樹結構的單元細分方案，對標記單元進行細分，返回步驟1。

圖2給出了本文特征值問題自適應分析流程。相較于原算法，主要區(qū)別以虛線流程框標注，即特征值誤差檢驗和超限單元標記。下文中，分別以縮寫APEF和APEV表示先后兩種不同的自適應方法(adaptive procedure guided by errors of eigen function / value)，后者即為本文所提出的算法。

圖2 特征值問題自適應分析流程Fig.2 Flow chart for adaptive analysis of eigenvalue problem

3 算法細節(jié)

3.1 有限元離散

采用有限元法對式(1)進行分析。在本文中，對于一維問題，采用m次多項式單元；對于二維問題，采用雙m次多項式單元，并結合基于四叉樹結構的單元細分法進行網(wǎng)格更新[23]，導出分層結構的網(wǎng)格(如圖3所示)。

圖3 二維問題分層網(wǎng)格劃分Fig.3 Hierarchical mesh refinement for 2D problem

(10)

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

確定單元及網(wǎng)格形式后，有限元求解歸結為求uh∈Sh滿足虛功方程

a(uh,vh)=λ·(uh,vh), ?vh∈Sh

(11)

經由標準的矩陣集成流程，無論一維或二維問題，原問題式(1)轉化為廣義矩陣特征值問題

KD=λMD

(12)

式中：D為振型向量；K和M分別為整體剛度矩陣和一致質量矩陣。

3.2 特征值計算

式(12)的廣義矩陣特征值問題，基于計數(shù)法和移位的逆冪(子空間)迭代，采用“劃界”和“求解”兩階段法求解，可保證不丟根、不落根，并已在APEF算法中得到了充分驗證；額外的精度導護措施，包括“μ檢驗”和“步長放大技術”，在本文中也同樣采納以確保解答的魯棒性。

特征值通過計算Rayleigh商獲得，可有效加速收斂過程。式(1)的各階次解將使得Rayleigh商，即式(13)取為駐值

(13)

(14)

3.3 特征值誤差估計

對式(13)取一階變分，有

(15)

取駐值條件，令一階變分為零，即得到式(3)。取二階變分，并利用式(3)，有

(16)

由級數(shù)展開，有

λ-λh=δλ+δ2λ+O(δ3λ)+…

(17)

由于δλ=0，故有

(18)

其中，δu=eh=u-uh，則有特征值的誤差估計

(19)

易知

a(eh,eh)≤C1h2m, (eh,eh)≤C2h2m+2

(20)

式中：h為單元尺寸參數(shù)；C1和C2為問題相關常數(shù)。故有Rayleigh商的誤差估計

|λ-λh|≤Ch2m

(21)

式中，C為問題相關常數(shù)?？芍邢拊獾腞ayleigh商誤差為2m階。

3.4 超收斂特征值計算

由式(21)，采用Rayleigh商計算特征值，其收斂階主要取決于位移和導數(shù)精度。若可提供收斂階更高的位移和導數(shù)，則理論上可獲得超收斂的特征值，也即

(22)

不失一般性地，本文以二階SL問題及彈性薄膜自由振動為模型問題進行討論，數(shù)學表達式分別為

(23a)

(23b)

式中：n為邊界外法向；ΓD為Dirichlet邊界；ΓN為Neumann邊界。

在算法設計時，有兩個方面需要著重說明。

3.4.1 單元積分策略

在先前基于EEP法的特征值問題自適應有限元分析研究中，采樣點主要用于單元上特征函數(shù)的誤差估計，因而按照圖4中采樣方式(1)均勻分布；不同于此，本文采樣點的另一作用在于作為積分點使用，以準確地計算式(22)，因而對于二維問題，應按照圖4中采樣方式(2)，也即相應數(shù)量Gauss點的位置分布。此外，對于本文采用的m次單元，通過EEP法獲得的超收斂解可能由m+2次的多項式表達。為保證積分準確性，保證在各方向上至少有m+2個采樣點。

圖4 單元采樣點分布Fig.4 Sampling points distribution over element domain

3.4.2 超收斂解的導數(shù)計算

在式(22)中，能量內積a(·,·)表達式中的積分項內一般包含位移的導數(shù)項。在本文中，采用對超收斂的EEP位移直接求導的方式進行計算。

對于SL問題，超收斂位移按式(24)計算，可直接求導

(24)

對于二維Laplace算子的特征值問題，兩個方向的超收斂位移導數(shù)需要區(qū)分計算。本文中超收斂位移的計算采用單元逐維精度恢復方案，表達式為

(25)

(26)

圖5 二維問題超收斂位移的導數(shù)計算Fig.5 Calculation of superconvergent derivatives for 2D problems

(27)

則對插值多項式PN-1(ξ)求導可得到ξ方向導數(shù)。

4 數(shù)值算例

本文算法已編制成Fortran 90程序。本章通過2個二階SL問題算例、2個彈性薄膜自由振動算例，驗證特征值超收斂計算及其驅動的自適應有限元分析的可靠性和有效性。

4.1 常截面桿軸向自由振動

本問題對應于式(23a)，參數(shù)p=r=1，q=0。對于該題，精確的特征值及特征函數(shù)表達為

λk=kπ2，uk=sin(kπx)

(28)

分別取低階和中高階的若干特征值，采用單元二分加密的方式，統(tǒng)計有限元解和超收斂解的收斂階。

如表1～表3所示，對于常系數(shù)問題，用有限元解uh計算和用EEP解u*計算特征值，后者收斂階可提升4階，即由2m階提升為2(m+2)階。

表1 常截面桿軸向自由振動第1階特征值收斂階結果(線性元)Tab.1 The convergence order of the 1st eigenvalue of the axial free vibration of a bar with constant cross section (Linear elements)

表2 常截面桿軸向自由振動第1階特征值收斂階結果(三次元)Tab.2 The convergence order of the 1st eigenvalue of the axial free vibration of a bar with constant cross section (Cubic elements)

表3 常截面桿軸向自由振動第5階特征值收斂階結果(二次元)Tab.3 The convergence order of the 5th eigenvalue of the axial free vibration of a bar with constant cross section (Quadratic elements)

4.2 局部震蕩問題

考慮式(29)參數(shù)定義的SL問題

p(x)=e10/cosh(10x),q(x)=0,
r(x)=100π2e10cosh(10x)/sinh210,
a=0,b=1,u(a)=0,u(b)=0

(29)

本例有精確解

(30)

該問題特點是在求解域(a,b)右端點附近的各階特征函數(shù)均呈現(xiàn)劇烈震蕩，求解有一定的困難。

首先驗證超收斂特征值的有效性。表4給出了該例第3階特征值的收斂階，對于該變系數(shù)問題，用EEP解求征值其精度可提升2階，即由2m階提升為2(m+1)階。

表4 震蕩問題第3階特征值收斂階結果(二次元)Tab.4 The convergence order of the 3rd eigenvalue of oscillation problem (Quadratic elements)

設誤差限為Tol=10-5，分別采用線性元及二次元，以APEV進行自適應分析。圖6～圖9對比了自適應及均勻網(wǎng)格加密的收斂過程，并繪制了單元尺寸分布?？梢钥吹?，基于本文自適應方案，收斂效率均有一定程度的提高，單元分布反映了問題的震蕩特性。

圖6 震蕩問題自適應及均勻加密過程(線性元)Fig.6 Adaptive and uniform mesh refinement for oscillation problem (Linear elements)

圖7 震蕩問題自適應網(wǎng)格單元尺寸分布(線性元)Fig.7 Element size distribution of adaptive mesh for oscillation problem (Linear elements)

為進一步說明APEV的特點，對本例設置更苛刻的誤差限Tol=10-8，作為對照，同時采用APEF進行計算。

采用三次元，從4個均分單元的初始網(wǎng)格開始，對第3階特征對進行自適應分析。圖10給出了兩種方法自適應過程的誤差下降圖?？梢钥吹剑捎谕瑫r控制特征值和特征函數(shù)誤差，APEF自適應方法存在較大的計算冗余度，算法停機時自由度數(shù)為947，特征值誤差為1.01×10-11，遠小于給定誤差限；而本文算法直接控制特征值誤差，算法停機時自由度數(shù)為308，特征值誤差為7.06×10-9，體現(xiàn)了較高的計算效率。

圖8 震蕩問題自適應及均勻加密過程(二次元)Fig.8 Adaptive and uniform mesh refinement for oscillation problem (Quadratic elements)

圖9 震蕩問題自適應網(wǎng)格單元尺寸分布(二次元)Fig.9 Element size distribution of adaptive mesh for oscillation problem (Quadratic elements)

圖10 震蕩問題新舊方法自適應過程Fig.10 Adaptive mesh refinement with current and previous methods for oscillation problem

4.3 固支扇形膜自由振動

考慮四邊固支的扇形膜(如圖11所示)，其模態(tài)對應于環(huán)形彈性膜的反對稱模態(tài)。本例中存在圓弧幾何邊界，具有一定的特殊性，采用精確幾何單元捕捉該特點。

圖11 扇形膜及三次幾何精確單元Fig.11 Fan-shaped membrane and cubic geometrically exact element

本例存在精確解。對于半徑為a≤r≤b的環(huán)形膜，其解析解為

(31)

式中，Jm和Ym分別對應于第一類和第二類的m階Bessel函數(shù)。解析的特征值表達式為

(32)

式中，km,n為式(33)的第n階根

(33)

對于本例，采用四次元，以單個單元網(wǎng)格為初始狀態(tài)，進行前8階特征值的自適應有限元分析，誤差限設為Tol=10-4，結果匯總于表5。誤差比顯示，各階特征值有限元解均滿足誤差要求。

表5 固支扇形膜自由振動特征值結果(APEV)Tab.5 Eigenvalues for free vibration of clamped fan-shaped membrane (APEV)

對本例，也采用APEF自適應策略進行分析，給定與以上相同的初始網(wǎng)格、單元及誤差限。表6給出了計算結果。與APEV相比，APEF的冗余度較高，以第1階為例，盡管使用了高達4 081個自由度，但其特征值結果反而差于APEV計算結果。

表6 固支扇形膜自由振動特征值結果(APEF)Tab.6 Eigenvalue for free vibration of clamped fan-shaped membrane (APEF)

圖12和圖13分別給出了第7階的自適應最終網(wǎng)格及有限元解模態(tài)。

圖12 固支扇形膜第7階最終網(wǎng)格Fig.12 7th final mesh of clamped fan-shaped membrane

圖13 固支扇形膜第7階模態(tài)Fig.13 7th mode for clamped fan-shaped membrane

為更充分說明自適應網(wǎng)格劃分的有效性，圖14給出了均勻網(wǎng)格下第7階模態(tài)的有限元解誤差分布，與圖12相一致地，自適應網(wǎng)格細密處誤差也較大，體現(xiàn)了自適應單元加密的合理性。

以第4階為例，圖15對自適應過程中特征值的收斂性進行了分析。圖15中，F(xiàn)EM和EEP分別表示在相同網(wǎng)格下經由有限元法和EEP后處理得到的特征值。在各個迭代步驟，基于超收斂特征值的誤差估計均逼近真實誤差，體現(xiàn)了其有效性；在停機時，誤差小于誤差限Tol，體現(xiàn)了APEV過程的有效性。

圖14 固支扇形膜均勻網(wǎng)格下第7階模態(tài)誤差Fig.14 Errors of 7th mode for clamped fan-shaped membrane on uniform mesh

圖15 固支扇形膜第4階特征值收斂分析Fig.15 Convergence analysis of the 4th order eigenvalues of clamped fan-shaped membrane

4.4 固支L形膜自由振動

考慮固支L形膜自由振動問題(如圖16所示)，其部分階模態(tài)位移導數(shù)在凹角處存在奇異性，采用常規(guī)有限元法具有較高的求解難度。采用三次元，以三個單元的網(wǎng)格作為初始狀態(tài)，誤差限設為Tol=10-3。

圖16 固支L形膜計算模型及初始網(wǎng)格Fig.16 Computational model and initial mesh for clamped L-shaped membrane

圖17給出了第1、第3、第5階的自適應最終網(wǎng)格及相應階次模態(tài)，本文算法自動識別了相應階次的模態(tài)特性。圖18及圖19給出了第1階和第5階特征值的收斂過程。對于本例，由于應力奇異點的存在，EEP法后處理特征值精度在自適應前期并不優(yōu)于有限元解，但局部誤差估計仍較為有效地體現(xiàn)了誤差分布特性，體現(xiàn)在網(wǎng)格可自動識別問題難點，向奇異點附近加密，并以與誤差限Tol相近的真實誤差停機。

圖17 固支L形膜計算模型及初始網(wǎng)格Fig.17 Meshes and modes for L-shaped membrane

圖18 固支L形膜第1階特征值收斂分析Fig.18 Convergence analysis of the 1st order eigenvalues of clamped L-shaped membrane

圖19 固支L形膜第5階特征值收斂分析Fig.19 Convergence analysis of the 5th order eigenvalues of clamped L-shaped membrane

5 結 論

本文建立了以頻率(特征值)誤差控制為目標的自適應有限元分析策略。以特征值誤差估計為核心，本文首先分析了Rayleigh商的收斂階，基于此提出了超收斂特征值的計算方案，其與原有限元解之差即可作為特征值的誤差估計。以特征值誤差估計作為整體停機標準，以特征函數(shù)誤差估計衡量局部單元誤差，本文通過兩次網(wǎng)格迭代間的誤差關系，將整體和局部的關系聯(lián)系起來，最終形成了完整自適應方案。

自適應分析的主要目標是采用“盡可能少”的自由度數(shù)獲得“盡可能高”的計算精度，其實現(xiàn)依賴于網(wǎng)格單元的合理分布，進一步依賴于對有限元解的有效誤差估計。本文數(shù)值結果顯示，對于光滑問題，相較于有限元解，EEP后處理特征值具有超收斂性，精度至少高2階；對于奇異問題，特征值盡管不具備超收斂性，APEV自適應分析過程仍具有有效性，顯著提升了有限元分析效率。此外，本文方法并不局限于作為模型的SL問題和彈性薄膜自由振動問題，在滿足分層結構網(wǎng)格的條件下，該方法可一般性地推廣至包括板、殼、三維彈性體在內的各類結構自由振動乃至彈性穩(wěn)定問題。