孫浩涵，袁 駟

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

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

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

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

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

1 背景知識

1.1 研究對象

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

Lu=λru

(1)

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

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

(2)

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

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

(3)

1.2 單元能量投影法

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

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

(4)

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

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

2 總體策略

2.1 自適應(yīng)分析目標(biāo)

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

(5)

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

(6)

2.2 整體和局部誤差

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

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

(7)

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

(8)

誤差大于該值的單元均應(yīng)細(xì)分。在實(shí)際算法中，為保證魯棒性，誤差限總?cè)闅v史路徑的最小值，即

(9)

2.3 整體算法

本文自適應(yīng)有限元分析總體策略可總結(jié)為：

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

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

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

步驟4單元細(xì)分。采用基于四叉樹結(jié)構(gòu)的單元細(xì)分方案，對標(biāo)記單元進(jìn)行細(xì)分，返回步驟1。

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

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

3 算法細(xì)節(jié)

3.1 有限元離散

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

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

(10)

式中：Ni和Nj為m次多項(xiàng)式形函數(shù)；dij為相應(yīng)的結(jié)點(diǎn)位移；單元檢驗(yàn)函數(shù)vh采用相同的插值形式。

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

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

(11)

經(jīng)由標(biāo)準(zhǔn)的矩陣集成流程，無論一維或二維問題，原問題式(1)轉(zhuǎn)化為廣義矩陣特征值問題

KD=λMD

(12)

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

3.2 特征值計(jì)算

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

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

(13)

(14)

3.3 特征值誤差估計(jì)

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

(15)

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

(16)

由級數(shù)展開，有

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

(17)

由于δλ=0，故有

(18)

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

(19)

易知

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

(20)

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

|λ-λh|≤Ch2m

(21)

式中，C為問題相關(guān)常數(shù)。可知有限元解的Rayleigh商誤差為2m階。

3.4 超收斂特征值計(jì)算

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

(22)

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

(23a)

(23b)

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

在算法設(shè)計(jì)時，有兩個方面需要著重說明。

3.4.1 單元積分策略

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

圖4 單元采樣點(diǎn)分布Fig.4 Sampling points distribution over element domain

3.4.2 超收斂解的導(dǎo)數(shù)計(jì)算

在式(22)中，能量內(nèi)積a(·,·)表達(dá)式中的積分項(xiàng)內(nèi)一般包含位移的導(dǎo)數(shù)項(xiàng)。在本文中，采用對超收斂的EEP位移直接求導(dǎo)的方式進(jìn)行計(jì)算。

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

(24)

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

(25)

(26)

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

(27)

則對插值多項(xiàng)式PN-1(ξ)求導(dǎo)可得到ξ方向?qū)?shù)。

4 數(shù)值算例

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

4.1 常截面桿軸向自由振動

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

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

(28)

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

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

表1 常截面桿軸向自由振動第1階特征值收斂階結(jié)果(線性元)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階特征值收斂階結(jié)果(三次元)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階特征值收斂階結(jié)果(二次元)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)

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

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

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

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

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

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

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

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

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

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

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

4.3 固支扇形膜自由振動

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

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

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

(31)

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

(32)

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

(33)

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

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

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

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

圖12和圖13分別給出了第7階的自適應(yīng)最終網(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

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

以第4階為例，圖15對自適應(yīng)過程中特征值的收斂性進(jìn)行了分析。圖15中，F(xiàn)EM和EEP分別表示在相同網(wǎng)格下經(jīng)由有限元法和EEP后處理得到的特征值。在各個迭代步驟，基于超收斂特征值的誤差估計(jì)均逼近真實(shí)誤差，體現(xiàn)了其有效性；在停機(jī)時，誤差小于誤差限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)位移導(dǎo)數(shù)在凹角處存在奇異性，采用常規(guī)有限元法具有較高的求解難度。采用三次元，以三個單元的網(wǎng)格作為初始狀態(tài)，誤差限設(shè)為Tol=10-3。

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

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

圖17 固支L形膜計(jì)算模型及初始網(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 結(jié) 論

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

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