平面曲梁面外自由振動(dòng)有限元分析的p 型超收斂算法

葉康生，梁 童

(清華大學(xué)土木工程系，土木工程安全與耐久教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084)

平面曲梁結(jié)構(gòu)在現(xiàn)代土木、航天、機(jī)械等領(lǐng)域具有廣泛的應(yīng)用，其自由振動(dòng)是結(jié)構(gòu)分析的重要內(nèi)容。該問(wèn)題目前的分析方法主要有：有限元法[1?3]、有限差分法[4]、微分容積法[5?6]、偽譜法[7]、動(dòng)力剛度法[8?10]等，其中有限元法應(yīng)用最為廣泛。

有限元法求解自由振動(dòng)問(wèn)題時(shí)，計(jì)算精度由單元次數(shù)和網(wǎng)格劃分決定。網(wǎng)格越密，單元次數(shù)越高，計(jì)算精度越好，但計(jì)算量亦隨之大幅攀升，高階頻率和振型尤甚[1]。因此，研究如何有效提高有限元求解自由振動(dòng)問(wèn)題的精度和效率具有十分重要的意義，超收斂計(jì)算則是解決這一難題的有效途徑。

目前結(jié)構(gòu)自由振動(dòng)問(wèn)題的各類(lèi)有限元超收斂求解方法，基本都是先對(duì)振型進(jìn)行超收斂修復(fù)，再將修復(fù)后的振型代入Rayleigh 商得到頻率的超收斂解，區(qū)別主要在于對(duì)振型的超收斂修復(fù)做法不同。這些方法中最典型的是文獻(xiàn)[11]提出的SPRD法，該法利用振型結(jié)點(diǎn)位移的超收斂性，在多個(gè)單元組成的小片上通過(guò)對(duì)振型用更高次多項(xiàng)式對(duì)結(jié)點(diǎn)位移進(jìn)行多項(xiàng)式擬合獲得振型的超收斂解。文獻(xiàn)[12]的做法與文獻(xiàn)[11]做法類(lèi)似，亦是基于振型結(jié)點(diǎn)位移的超收斂性，在幾個(gè)單元構(gòu)成的小片上通過(guò)對(duì)振型用更高次多項(xiàng)式對(duì)結(jié)點(diǎn)位移作多項(xiàng)式插值獲得振型的超收斂解；文獻(xiàn)[13?14]僅修復(fù)振型所對(duì)應(yīng)的應(yīng)力，并代入Rayleigh 商的應(yīng)變能部分，通過(guò)應(yīng)變能計(jì)算精度的提高獲得頻率的超收斂解；文獻(xiàn)[15?16]利用等幾何分析法，通過(guò)構(gòu)造高階質(zhì)量矩陣，獲得頻率的高精度解；另外還有文獻(xiàn)[17]基于外推的方法和文獻(xiàn)[18]基于投影的方法。

文獻(xiàn)[19]對(duì)自由振動(dòng)問(wèn)題提出一種p型超收斂算法，結(jié)果顯示，該法簡(jiǎn)單、高效，能顯著提升頻率和振型的精度和收斂階。文獻(xiàn)[20]將該法應(yīng)用于平面曲梁的靜力分析中，文獻(xiàn)[21]將該法延拓至平面曲梁的面內(nèi)自由振動(dòng)分析中，文獻(xiàn)[22]將該法延拓至歐拉梁的彈性穩(wěn)定分析中，文獻(xiàn)[23]將該法延拓至非線性常微分方程邊值問(wèn)題的分析中。本文的工作是文獻(xiàn)[19?23]工作的延續(xù)，是文獻(xiàn)[21]工作的后半部分，進(jìn)一步將p型超收斂算法應(yīng)用于平面曲梁的面外自由振動(dòng)分析中。數(shù)值算例表明，該法仍延續(xù)了求解前述問(wèn)題時(shí)的優(yōu)秀表現(xiàn)，算法穩(wěn)定、高效，能顯著提高解答的精度、質(zhì)量和收斂階，再次展現(xiàn)了p型超收斂算法在求解自由振動(dòng)問(wèn)題時(shí)的優(yōu)良特性。

1 平面曲梁自由振動(dòng)分析模型

平面曲梁的軸線為平面曲線，所有截面關(guān)于軸線所在平面對(duì)稱(chēng)。如圖1 所示，建立曲梁分析的xyz坐標(biāo)系，x、y分別沿軸線切向、法向，z垂直于軸線所在平面。記s為軸線弧長(zhǎng)坐標(biāo)。曲梁面外自由振動(dòng)時(shí)，分別記面外位移、y軸轉(zhuǎn)角位移和扭轉(zhuǎn)角位移為w(s) 、 θ(s) 、 ?(s) ，對(duì)應(yīng)內(nèi)力Q、M、T如圖1 所示。

圖1 平面曲梁Fig.1 Planar curved beam

面外自由振動(dòng)時(shí)，曲梁的幾何方程、物理方程和平衡方程分別為[24]：

式中： ()′=d()/ds；E為彈性模量；G為剪切模量； ρ為材料密度；A為截面面積；Iy為截面關(guān)于y軸的慣性矩；Ip為截面極慣性矩；k為截面形狀系數(shù)；J為截面扭轉(zhuǎn)常數(shù)；R為軸線曲率半徑，當(dāng)軸線曲率中心位于y軸正向時(shí)R為正，位于y軸負(fù)向時(shí)R為負(fù)。

將物理方程和幾何方程代入平衡方程，得到用位移表示的平面曲梁面外自由振動(dòng)的控制微分方程為：

式中：l為曲梁的軸線長(zhǎng)度；L、m為相應(yīng)的微分算子矩陣，其具體表達(dá)式可由式(4)導(dǎo)出，此處不再詳細(xì)給出。該問(wèn)題理論上有無(wú)窮多組解答(λn,dn(s)),n=1,2,··· ，特征值λn按從小到大排序：

用有限元法求解該問(wèn)題時(shí)，先對(duì)結(jié)構(gòu)進(jìn)行有限元離散，將結(jié)構(gòu)劃分為ne個(gè)單元，如圖2 所示。

記網(wǎng)格尺寸：

圖2 有限元網(wǎng)格劃分Fig.2 Finite element mesh

式中，he為單元 (e) 的軸線長(zhǎng)度，he=se+1?se。雖然理論上d={w,θ,?}T中三個(gè)位移分量的形函數(shù)階數(shù)可取為不同，但數(shù)值結(jié)果表明，有限元解的收斂階取決于這三個(gè)位移分量形函數(shù)階數(shù)中的最低值，為使有限元自由度的求解效率得到充分發(fā)揮，應(yīng)讓三個(gè)位移分量的形函數(shù)階數(shù)取為相同，故本文有限元計(jì)算和超收斂計(jì)算中三個(gè)位移分量的形函數(shù)階數(shù)均取為相同。在任意單元 (e)內(nèi)，采用m次Hierarchical 形函數(shù)將解答近似為：

式中： ξ=(s?se)/he，為單元 (e)上的無(wú)量綱坐標(biāo)；we、 θe、 ?e為第e個(gè)結(jié)點(diǎn)上的結(jié)點(diǎn)位移；ai、bi、ci(i=1,2,···,m?1) 分別為w、 θ 、 ?的單元內(nèi)部自由度。將式(9)寫(xiě)成矩陣形式：

式中， ?e為單元 (e)的自由度向量，具體為：

Nw、Nθ、N?分別為w、 θ 、 ?的形函數(shù)矩陣，其具體表達(dá)式可由式(9)和式(11)導(dǎo)出。

由變分原理[1]，有限元法將原問(wèn)題離散為如下矩陣廣義特征值問(wèn)題：

式中：Kh和Mh分別為該網(wǎng)格下的整體剛度矩陣和整體質(zhì)量矩陣，分別由單元?jiǎng)偠染仃嚭蛦卧|(zhì)量矩陣集成獲得； ?h為該網(wǎng)格下的結(jié)構(gòu)整體自由度向量。單元?jiǎng)偠染仃嚭蛦卧|(zhì)量矩陣具體計(jì)算如下：

求解上述矩陣廣義特征值問(wèn)題即可獲得該網(wǎng)格下的有限元解 (λhn,dnh(s))，n=1,2,···。

對(duì)該向量型常微分方程特征值問(wèn)題，文獻(xiàn)[1]估計(jì)有限元解中頻率和振型分別具有如下收斂階：

式中，m為單元多項(xiàng)式次數(shù)。

文獻(xiàn)[25]對(duì)該向量型問(wèn)題對(duì)應(yīng)的靜力問(wèn)題Ld=f只有一個(gè)分量的最簡(jiǎn)單情形，證明當(dāng)采用m次元求解時(shí)，結(jié)點(diǎn)位移具有h2m階的超收斂性。數(shù)值算例表明，該結(jié)論對(duì)特征值問(wèn)題同樣成立，在自由振動(dòng)問(wèn)題中，振型結(jié)點(diǎn)位移同樣具有h2m階的超收斂性，即：

振型結(jié)點(diǎn)位移的超收斂性是本文算法的基礎(chǔ)。

2 超收斂求解思路

由式(5)可知，在單元 (e) 上，第n階振型dn(s)滿(mǎn)足如下常微分方程邊值問(wèn)題：

且當(dāng)網(wǎng)格足夠密時(shí)，為該邊值問(wèn)題的唯一解。

由式(15)和式(17)知，頻率和振型結(jié)點(diǎn)位移均具有h2m階的超收斂性。故在網(wǎng)格加密過(guò)程中，這些值將比內(nèi)點(diǎn)位移更快地向精確解收斂，當(dāng)網(wǎng)格足夠密時(shí)，有：

此時(shí)，單元 (e)上的振型將近似于下列線性邊值問(wèn)題的解：

3 誤差估計(jì)

4 數(shù)值算例

本文算法已編成Fortran 程序。下面通過(guò)3 個(gè)數(shù)值算例來(lái)展示本文算法的計(jì)算精度和修復(fù)效果。為簡(jiǎn)單起見(jiàn)，三個(gè)算例均采用無(wú)量綱計(jì)算。為考察振型收斂階，本文選取振型誤差向量模為長(zhǎng)度的最大模：

平面曲梁面外自由振動(dòng)問(wèn)題只有常截面圓弧曲梁具有解析解。為展現(xiàn)本文算法對(duì)收斂階的提升，例1 和例2 選取了不同邊界條件的兩個(gè)常截面圓弧曲梁；為展現(xiàn)本文算法的修復(fù)效果，例3 與已有文獻(xiàn)的計(jì)算結(jié)果進(jìn)行了對(duì)比。

例1. 常截面圓弧曲梁1

圖3 例1 常截面圓弧曲梁Fig.3 Uniform arch beam in Example 1

可求出該圓弧曲梁頻率和振型的精確解，根據(jù)文獻(xiàn)[26]，設(shè)其振型為：

式中，n≥0且為整數(shù)。將式(42)代入式(4)，化簡(jiǎn)可得關(guān)于 {w0,θ0,?0}T的齊次方程組。由于振型存在非零解，故該齊次方程組的系數(shù)矩陣奇異，其行列式值為零，從而可推出如下的頻率方程：

由該頻率方程可求出每個(gè)n值對(duì)應(yīng)的三個(gè)頻率值。本例第一階、第二階頻率均為0，故對(duì)第三階頻率求解。第三階頻率對(duì)應(yīng)n=2，頻率方程如式(44)所示，求解該方程可得如式(45)所示頻率。將求出的頻率代入 {w0,θ0,?0}T的齊次方程組，可求得對(duì)應(yīng)的振型，如式(46)所示，該階振型取w(α=0)=1進(jìn)行歸一化。

首先考察頻率、振型和振型結(jié)點(diǎn)位移有限元解的超收斂性。對(duì)該階振型分別采用二次元(m=2) 、三次元 (m=3) 和四次元 (m=4)進(jìn)行有限元求解，計(jì)算頻率誤差 |???h|， 振型誤差||d?dh||和跨中截面結(jié)點(diǎn)C處( α=π/2 )位移誤差 ||dc?dch||并進(jìn)行收斂階的統(tǒng)計(jì)，結(jié)果如圖4 和表1 所示。結(jié)果表明頻率和振型結(jié)點(diǎn)位移具有 2m階的超收斂性，而振型的收斂階僅為m+1階。

圖4 例1 第三階頻率、振型和結(jié)點(diǎn)位移的收斂階Fig.4 Rate of convergence of 3rd frequence, mode and its nodal displacement in Example 1

為考察p型超收斂算法的計(jì)算精度和修復(fù)效果，將曲梁均勻劃分為5 個(gè)單元，先用一次元(m=1) 進(jìn)行有限元求解，再用二次元 (=2)進(jìn)行超收斂修復(fù)。得到第三階頻率的有限元解為?h3=8.398 (相對(duì)誤差為300.5%)，超收斂解為??3=2.664(相對(duì)誤差為27.0%)。圖5 為第三階振型的有限元解、超收斂解和精確解的對(duì)比。可見(jiàn)，本文算法可以顯著提升頻率和振型的計(jì)算精度。

為考察p型超收斂算法的收斂階，對(duì)第三階頻率和振型，先用二次元 (m=2)進(jìn)行有限元計(jì)算，再用三次元 (=3) 、四次元 (=4)和五次元(=5)進(jìn)行超收斂修復(fù)并統(tǒng)計(jì)收斂階，計(jì)算結(jié)果如圖6 和表2 所示。易見(jiàn)，超收斂修復(fù)后，精度和收斂階與有限元解相比均有顯著提升。同時(shí)，由表2 可以看出，當(dāng)從四次元 (=4)升至五次元(=5)時(shí)，頻率和振型的收斂階未隨之提升。理論分析和數(shù)值算例均表明，當(dāng)>2m時(shí)，超收斂

解的收斂階與=2m時(shí)相同，不會(huì)隨的進(jìn)一步提高而增加。

表1 例1 第三階頻率、振型和結(jié)點(diǎn)位移的收斂階Table 1 Rate of convergence of 3rd frequency, mode and its nodal displacement in Example 1

圖5 例1 第三階振型的有限元解與超收斂解比較Fig.5 Comparison of FE solution with recovered solution on 3rd mode in Example 1

圖6 例1 第三階特征對(duì)有限元解與超收斂解收斂階Fig.6 Rate of convergence of FE solution and recovered solutions on 3rd eigenpair in Example 1

為對(duì)比本文算法與SPRD 法[11]的計(jì)算效果，對(duì)第三階頻率和振型的二次元解用SPRD 法進(jìn)行超收斂修復(fù)，用相鄰四個(gè)結(jié)點(diǎn)(單元兩端結(jié)點(diǎn)外加離單元最近的兩個(gè)結(jié)點(diǎn))上的位移解作三次多項(xiàng)式插值得單元上振型的超收斂解，將振型超收斂解代入Rayleigh 商得頻率超收斂解。頻率和振型的收斂階統(tǒng)計(jì)結(jié)果如表3 所示?？梢?jiàn)，SPRD 法三次元修復(fù)解與本文算法三次元修復(fù)解的收斂階相同，只是同樣網(wǎng)格下SPRD 法的解答精度比本文算法的解答精度稍低些。

例2. 常截面圓弧曲梁2

本例仍取軸線為半個(gè)圓周的常截面圓弧曲梁，曲梁參數(shù)同例1，兩端支座形式改為面外簡(jiǎn)支，如圖7 所示。

其振型可設(shè)為：

表2 例1 第三階特征對(duì)有限元解與超收斂解收斂階Table 2 Rate of convergence of FE solution and recovered solutions on 3rd eigenpair in Example 1

式中，n≥0且為整數(shù)。

將曲梁均勻劃分為32 個(gè)單元 (ne=32)，先用一次元 (m=1)進(jìn)行有限元求解，再用二次元(=2)進(jìn)行超收斂修復(fù)。前12 階頻率和振型的有限元解與超收斂解結(jié)果如表4 所示。易見(jiàn)，對(duì)一次元，雖然振型結(jié)點(diǎn)位移不具有超收斂性，本文方法仍能顯著提高解答(尤其是頻率)的精度。

表3 例1 第三階特征對(duì)SPRD 法收斂階Table 3 Rate of convergence of recovered solutions on 3rd eigenpair in Example 1 with SPRD method

為考察本文算法對(duì)高階頻率和振型的計(jì)算效果，對(duì)該曲梁第26 階頻率和振型進(jìn)行分析，此階對(duì)應(yīng)n=14。按前述方法計(jì)算該階頻率和振型的精確解，精確解為：

圖7 例2 常截面圓弧曲梁Fig.7 Uniform arch beam in Example 2

先用三次元 (m=3)進(jìn)行有限元求解，再依次用四次元 (=4) 、五次元 (=5) 和六次元(=6)進(jìn)行超收斂修復(fù)，結(jié)果如圖8 和表5 所示。結(jié)果表明p型超收斂算法對(duì)曲梁面外自由振動(dòng)的高階頻率和振型仍有顯著的修復(fù)效果。

表4 例2 前12 階特征對(duì)有限元解與超收斂解誤差Table 4 Errors of FE solution and recovered solutions on first 12 eigenpairs in Example 2

例3. 一般軸線曲梁

文獻(xiàn)[27?28]研究了圖9 所示的拋物線、三角正弦線和橢圓線三類(lèi)軸線曲梁的面外自由振動(dòng)問(wèn)題，三類(lèi)曲梁軸線方程分別如下：

圖8 例2 第26 階特征對(duì)有限元解與超收斂解收斂階Fig.8 Rate of convergence of FE solution and recovered solutions on 26th eigenpair in Example 2

將三類(lèi)軸線曲梁沿軸線水平投影方向均分為32 個(gè)單元 (ne=32) ，先采用一次元 (m=1)進(jìn)行有限元計(jì)算，后采用二次元 (=2)進(jìn)行超收斂修復(fù)和先采用二次元 (m=2)進(jìn)行有限元計(jì)算，后采用三次元 (=3)進(jìn)行超收斂修復(fù)，計(jì)算前3 階頻率如表6 所示。與文獻(xiàn)結(jié)果相比，頻率在一次元有限元計(jì)算后仍誤差較大，但采用二次元進(jìn)行超收斂修復(fù)后，誤差顯著降低。而用二次元有限元計(jì)算后采用三次元進(jìn)行超收斂修復(fù)，得到的解答與文獻(xiàn)結(jié)果已十分接近。表7 給出了上述網(wǎng)格下三類(lèi)軸線曲梁前20 階振型有限元計(jì)算和超收斂計(jì)算的總耗時(shí)?？梢?jiàn)，p型超收斂計(jì)算耗時(shí)很少，但解答精度提升非常顯著。在一次元有限元計(jì)算后用二次元進(jìn)行超收斂修復(fù)，所得結(jié)果與直接用二次元進(jìn)行有限元計(jì)算的精度已很相近，但計(jì)算總耗時(shí)不足后者的1/2，充分展現(xiàn)出p型超收斂算法提升精度的高效性。

綜合數(shù)值算例，本文超收斂解的收斂規(guī)律如表8 所示。

表5 例2 第26 階特征對(duì)有限元解與超收斂解收斂階Table 5 Rate of convergence of FE solution and recovered solutions on 26th eigenpair in Example 2

表6 例3 三類(lèi)軸線曲梁的無(wú)量綱頻率Table 6 Dimensionless frequencies for three curved beams of Example 3

圖9 三類(lèi)軸線曲梁Fig.9 Three types of curved beams

表7 例3 前20 階振型計(jì)算總耗時(shí) /sTable 7 Total computation time for first 20 modes of Example 3

表8 超收斂解的收斂階Table 8 Convergence order of recovered solutions

5 結(jié)論

本文的p型超收斂算法在有限元解的基礎(chǔ)上，充分利用頻率和振型結(jié)點(diǎn)位移所具有的超收斂特性，通過(guò)逐單元求解振型近似滿(mǎn)足的線性邊值問(wèn)題獲得頻率和振型的更高精度解，方法思路巧妙，穩(wěn)定高效，能顯著提高解答精度和收斂階，可進(jìn)一步推廣應(yīng)用到包括薄壁結(jié)構(gòu)等空間結(jié)構(gòu)自由振動(dòng)的分析中。