王永亮

(1.中國(guó)礦業(yè)大學(xué)(北京)力學(xué)與建筑工程學(xué)院，北京100083；2.中國(guó)礦業(yè)大學(xué)(北京)煤炭資源與安全開(kāi)采國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京100083)

梁構(gòu)件在實(shí)際工程中廣泛應(yīng)用，如結(jié)構(gòu)工程中的曲線梁橋[1]、航空航天工程中的機(jī)翼翼梁[2]、采礦工程中的礦柱Euler–Bernoulli梁模型[3]等。梁的自由振動(dòng)是強(qiáng)迫振動(dòng)研究的基礎(chǔ)，自振頻率可為避免結(jié)構(gòu)共振提供依據(jù)[4]，利用振型疊加可分析地震載荷作用下結(jié)構(gòu)的變形量[5]。此外，利用自振頻率和振型可以進(jìn)行含損傷結(jié)構(gòu)的損傷識(shí)別[6?7]，損傷個(gè)數(shù)、位置和程度的準(zhǔn)確識(shí)別依賴(lài)于高精度的頻率和振型解答[8?9]。本文針對(duì)典型的平面變截面變曲率梁構(gòu)件自由頻率和振型求解開(kāi)展研究，該類(lèi)曲梁的軸線在同一平面內(nèi)、且曲梁的任意位置橫截面均關(guān)于上述平面對(duì)稱(chēng)。這類(lèi)曲梁在面內(nèi)振動(dòng)和面外振動(dòng)相互解耦，需要分別求解[10]，如何精確、有效地求解變截面變曲率梁的連續(xù)階頻率和振型成為需求，也是本文的研究目標(biāo)。

有限元法是求解曲梁自由振動(dòng)頻率和振型近似解的重要方法，被應(yīng)用于分析不同橫截面[11]、支撐類(lèi)型[12]、曲線線型形式[13]等復(fù)雜梁構(gòu)件的自由振動(dòng)問(wèn)題。以控制和提高有限元解的精度為目標(biāo)，自適應(yīng)有限元法被應(yīng)用于Euler-Bernoulli梁構(gòu)件[6]、梁系結(jié)構(gòu)[14?15]的頻率和振型的求解，各類(lèi)新型網(wǎng)格自適應(yīng)算法[16? 17]為提供滿(mǎn)足預(yù)設(shè)誤差限的高精度解答提供了高性能計(jì)算方案。本文將建立變截面變曲率梁面內(nèi)和面外自由振動(dòng)振型的有限元超收斂拼片恢復(fù)解的計(jì)算方法，并基于振型超收斂解構(gòu)建網(wǎng)格自適應(yīng)分析策略，可獲得優(yōu)化的網(wǎng)格和滿(mǎn)足預(yù)設(shè)誤差限的連續(xù)階頻率和振型。數(shù)值算例檢驗(yàn)了本文方法分析不同曲線形態(tài)、多類(lèi)邊界條件、變曲率、變化截面形式的曲梁面內(nèi)和面外自由振動(dòng)連續(xù)階頻率和振型的求解效力，表明本文方法分析結(jié)果的精確性和有效性。

1 曲梁面內(nèi)和面外自由振動(dòng)

考慮圖1所示的平面變截面變曲率梁，曲梁的中性軸坐標(biāo)為s，坐標(biāo)系為xy z，其中x、y為曲梁平面內(nèi)坐標(biāo)，x沿軸線切向，y沿軸線法向，z垂直于軸線所在平面。面內(nèi)振動(dòng)的位移為：沿x軸位移振幅u、沿y軸位移振幅v和繞z軸的轉(zhuǎn)角振幅ψz；面外振動(dòng)的位移為：沿z軸位移振幅w、繞y軸的轉(zhuǎn)角振幅 ψy和繞x軸的扭轉(zhuǎn)角振幅φx。記曲梁曲率半徑為R(s)，截面剪切剛度修正系數(shù)為κ，截面面積為A(s)，對(duì)y軸慣性矩為Iy(s)，對(duì)z軸慣性矩為I z(s)，極慣性矩為Ip(s)，扭轉(zhuǎn)常數(shù)為J(s)，長(zhǎng)度為l。記材料彈性模量為E，剪切模量為G，泊松比為ν，密度為 ρ。

圖1 平面變截面變曲率梁坐標(biāo)系和符號(hào)Fig.1 Coordinatesystemsand symbols of planar nonuniform and variablecurvature curved beam

本文研究的平面曲梁面內(nèi)自由振動(dòng)的微分控制方程[18]為：

式中：()′=d()/ds；ω為自振頻率。

面內(nèi)、面外自由振動(dòng)控制方程式(1)、式(2)均可記為如下矩陣形式的特征值方程：

式中：u為對(duì)應(yīng)的振型(位移)函數(shù)向量(面內(nèi)、面外振動(dòng)分別為{u,v,ψz}T、{w,ψy,φx}T)，ω、u分別對(duì)應(yīng)于特征值、特征向量，本文亦將(ω,u)合稱(chēng)為特征對(duì)；L、R是相應(yīng)的微分算子矩陣。

2 振型有限元解

對(duì)于求解特征值方程式(3)，在給定有限元網(wǎng)格π下，常規(guī)有限元建立如下線性矩陣特征值方程：

式中：D為振型向量的有限元解；K和M是靜力剛度矩陣和一致質(zhì)量矩陣。采用逆冪迭代法[20]求解，即得當(dāng)前網(wǎng)格下有限元解(ωh,u h)。

3 振型超收斂拼片恢復(fù)解

位移型有限元法以節(jié)點(diǎn)位移為基本未知量，有限元位移解在單元內(nèi)部具有hm+1階(m為當(dāng)前形函數(shù)階次)的超收斂性，在單元節(jié)點(diǎn)上具有h2m階的超收斂性，這些節(jié)點(diǎn)成為位移超收斂點(diǎn)。超收斂點(diǎn)上解答的誤差比其他區(qū)域解答更小，具有更快的收斂速度[20]。利用相鄰單元拼片，并將這些單元節(jié)點(diǎn)位移值進(jìn)行次高階形函數(shù)(高于當(dāng)前形函數(shù)階次m，即≥m+1)插值，即可獲得單元域內(nèi)位移場(chǎng)的超收斂解，形成有限元后處理超收斂拼片恢復(fù)方法[6,21?22]。在當(dāng)前網(wǎng)格下，位移超收斂解比常規(guī)有限元解具有更高收斂階。本文對(duì)于曲梁的面內(nèi)和面外自由振動(dòng)問(wèn)題，求得當(dāng)前網(wǎng)格下振型(位移)的有限元解后，利用有限元后處理超收斂拼片恢復(fù)方法，將單元e的相鄰單元進(jìn)行拼片，并進(jìn)行高階形函數(shù)插值處理，可以得到振型的超收斂解：

式中，a( )、b( )為應(yīng)變能、動(dòng)能內(nèi)積。

運(yùn)用Rayleigh 商計(jì)算得到的頻率超收斂解比振型超收斂解具有更高的收斂階[23]，因此，本文針對(duì)振型進(jìn)行誤差估計(jì)和控制，即可獲得高精度的振型解，并確保頻率解的準(zhǔn)確性。

4 誤差估計(jì)和網(wǎng)格細(xì)分加密

利用振型誤差估計(jì)，網(wǎng)格可以進(jìn)行優(yōu)化處理來(lái)降低和控制振型的誤差，達(dá)到預(yù)設(shè)的解答精度。本文方法采用式(14)對(duì)每個(gè)有限元單元e上的振型誤差進(jìn)行判斷，如果誤差控制式(14)不滿(mǎn)足，則表明該單元上振型解答的誤差過(guò)大，需要通過(guò)進(jìn)行網(wǎng)格優(yōu)化處理，本文采用單元均勻細(xì)分加密的方式來(lái)增加模型自由度、降低單元上解答的誤差[6,24]。當(dāng)前單元細(xì)分生成的新單元長(zhǎng)度與目前誤差和單元階次相關(guān)，即利用當(dāng)前誤差可以估計(jì)新單元長(zhǎng)度：

為適應(yīng)問(wèn)題的復(fù)雜性，往往在求解域上使用非均勻分布網(wǎng)格才能高效地獲得可靠解答。上述自適應(yīng)分析過(guò)程將單元均勻細(xì)分加密，如果每次細(xì)分出過(guò)多單元，則容易造成的最終單元的冗余，增大計(jì)算規(guī)模。本文方法對(duì)每次單元細(xì)分加密數(shù)目進(jìn)行控制，單元均勻細(xì)分的實(shí)際數(shù)目為：

5 自適應(yīng)分析策略

對(duì)曲梁自由振動(dòng)問(wèn)題，利用上述基于振型超收斂拼片恢復(fù)解的有限元解誤差估計(jì)、網(wǎng)格細(xì)分加密方法，可以形成有效的網(wǎng)格自適應(yīng)分析方法，最終可以得到優(yōu)化的網(wǎng)格，并獲得該優(yōu)化網(wǎng)格上滿(mǎn)足誤差限的各階高精度頻率和振型解答。對(duì)于振型，超收斂解u?較當(dāng)前網(wǎng)格而下有限元解u h具有更高的收斂階，即u?比u h更接近精確解u，因此采用u?代u對(duì)u h進(jìn)行誤差估計(jì)和控制，形成能量模形式的誤差估計(jì)；對(duì)于頻率，利用位移超收斂解計(jì)算Rayleigh 商[23]得到頻率的超收斂解。通過(guò)Rayleigh 商得出的頻率超收斂解比振型超收斂解具有更高的收斂階，本文通過(guò)估計(jì)和控制振型解答的誤差可保證頻率解答滿(mǎn)足誤差要求，如此可形成本文算法通過(guò)控制振型誤差的停機(jī)準(zhǔn)則?；谟邢拊獾恼`差估計(jì)進(jìn)行網(wǎng)格細(xì)分加密，即可求得優(yōu)化的網(wǎng)格和滿(mǎn)足誤差限的各階高精度頻率和振型解答，形成一套有限元自適應(yīng)求解方案：

1)有限元解：對(duì)于曲梁振動(dòng)問(wèn)題，在當(dāng)前網(wǎng)格π 上進(jìn)行常規(guī)有限元計(jì)算，得到當(dāng)前網(wǎng)格下特征對(duì)的有限元解(ωh,u h)；

2)振型超收斂解：利用有限元后處理超收斂拼片恢復(fù)方法，可得當(dāng)前網(wǎng)格下振型的超收斂解u?；同時(shí)，可得超收斂解與有限元解之間的誤差估計(jì)值ξ；

3)網(wǎng)格細(xì)分加密：對(duì)于 ξ不滿(mǎn)足式誤差控制式(14)的單元，用網(wǎng)格細(xì)分加密方法將其細(xì)分，得到新的有限元網(wǎng)格，返回步驟1)；如果所有單元均滿(mǎn)足誤差限，則網(wǎng)格無(wú)需再細(xì)分，求解過(guò)程結(jié)束。

6 數(shù)值算例

本文方法已經(jīng)編制Fortran 90程序，本節(jié)給出求解多種變截面變曲率梁面內(nèi)和面外自由振動(dòng)的數(shù)值算例，通過(guò)對(duì)比典型的拋物線曲梁、變截面變曲率梁、橢圓弧曲梁面內(nèi)自由振動(dòng)和拋物線曲梁、圓弧曲梁面外自由振動(dòng)結(jié)果來(lái)驗(yàn)證本文方法計(jì)算結(jié)果的精確性和適用性。本節(jié)所有算例均采用3次元，初始網(wǎng)格采用2 個(gè)單元，給定的初始誤差限為T(mén)ol=10?4，設(shè)定單元細(xì)分的極限數(shù)目為d=6。

例1.拋物線曲梁面內(nèi)自由振動(dòng)

考慮如圖2所示的拋物線曲梁，曲梁跨度L=5 m、高度h=4 m。

圖2 拋物線曲梁幾何模型Fig.2 Geometric model of parabolic curved beam

該曲線梁的拋物線方程為：

該曲線梁的基本參數(shù)如下：

表1 拋物線曲梁面內(nèi)自由振動(dòng)無(wú)量綱頻率值Table 1 Non-dimensional frequencies of in-plane vibration of parabolic curved beam

圖3 本文方法解答與密集網(wǎng)格高精度解的誤差Fig.3 Errorsof results of proposed method and highprecision solutionsunder dense mesh

圖4～圖6分別給出了本文方法求解兩端固支邊界、高跨比h/L=0.8情況下的前3階振型結(jié)果，并在橫坐標(biāo)軸上標(biāo)記出自適應(yīng)網(wǎng)格的最終分布情況。為方便直觀顯示和對(duì)比分析，圖中振型結(jié)果均進(jìn)行歸一化處理(令最大振型值為1)?？梢钥闯觯裥驮趦晒潭ǘ说奈灰凭鶠?值；本文方法求解各階振型均劃分出非均勻網(wǎng)格，且在振型變化平緩區(qū)域使用稀疏網(wǎng)格、在振型變化劇烈處采用了相對(duì)細(xì)密的網(wǎng)格，避免了全域使用一致細(xì)密網(wǎng)格的冗余性。

圖4 第1階振型(uh1,vh1,ψ hz1)Fig.4 First order vibration mode(uh1,vh1,ψ hz1)

圖5 第2階振型(uh2,vh2,ψhz2)Fig.5 Second order vibration mode (uh2,vh2,ψ hz2)

圖6 第3階振型Fig.6 Third order vibration mode (uh3,vh3,ψ hz3)

例2.變截面變曲率梁面內(nèi)自由振動(dòng)

考慮圖7所示的變截面變曲率梁，曲梁兩端為固支邊界條件、跨度為L(zhǎng)。

圖7 變截面變曲率梁幾何模型Fig.7 Geometric model of non-uniform and variable curvaturecurved beam

使用本文方法連續(xù)求解該曲梁面內(nèi)自由振動(dòng)的連續(xù)前4階特征對(duì)，計(jì)算頻率值列于表2。文獻(xiàn)[26]發(fā)展基于Wittrick-Williams算法的動(dòng)力剛度法求解該問(wèn)題，解答如表2所示，可以看出本文求解結(jié)果與該動(dòng)力剛度法結(jié)果吻合較好，檢驗(yàn)了本文方法對(duì)同時(shí)變化截面和曲率形式梁的面內(nèi)自由振動(dòng)求解的有效性。

表2 變截面變曲率梁面內(nèi)自由振動(dòng)頻率值Table 2 Frequencies of in-plane vibration of non-uniform and variable curvature curved beam

例3.橢圓弧曲梁面內(nèi)自由振動(dòng)

例4.拋物線曲梁面外自由振動(dòng)

考慮例1中圖2所示拋物線曲梁的面外自由振動(dòng)，拋物線方程同式(18)，該曲線梁的基本參數(shù)為：

圖8 橢圓弧曲梁幾何模型Fig.8 Geometric model of of elliptic curved beam

表3 橢圓弧曲梁面內(nèi)自由振動(dòng)無(wú)量綱頻率值Table 3 Non-dimensional frequenciesof in-plane vibration of elliptic curved beam

表4 拋物線曲梁面外自由振動(dòng)無(wú)量綱頻率值Table 4 Non-dimensional frequencies of out-of-plane vibration of parabolic curved beam

表5 圓弧曲梁面外自由振動(dòng)無(wú)量綱頻率值Table5 Non-dimensional frequenciesof out-of-plane vibration of circular curved beam

例5.圓弧曲梁面外自由振動(dòng)

考慮圖10所示兩端固支圓弧曲梁，圓弧角度為 θ0，梁橫截面為圓形。

圖10 圓弧曲梁幾何模型Fig.10 Geometric model of circular curved beam

該曲線梁的基本參數(shù)為：

7 結(jié)論

本文研究工作的主要結(jié)論和計(jì)算方法潛力可總結(jié)如下：

(1)本文提出變截面變曲率梁面內(nèi)和面外自由振動(dòng)振型的有限元后處理超收斂拼片恢復(fù)方法，建立各階振型的超收斂解，基于振型超收斂解進(jìn)行網(wǎng)格自適應(yīng)加密，求解得到優(yōu)化的網(wǎng)格與滿(mǎn)足用戶(hù)給定誤差限的頻率和振型。

(2)本文算法具有通用性，分析了不同曲線形態(tài)、多類(lèi)邊界條件、變截面、變曲率形式的曲梁振動(dòng)問(wèn)題，可推廣到一般形式結(jié)構(gòu)和固體變形場(chǎng)有限元解的超收斂計(jì)算及自適應(yīng)分析。

(3)本文方法具有應(yīng)用分析含裂紋損傷曲線形桿件振動(dòng)和穩(wěn)定等工程實(shí)際問(wèn)題的潛力，獲得損傷缺陷下的高精度振型和失穩(wěn)模態(tài)，為分析裂紋損傷對(duì)曲線形桿系結(jié)構(gòu)的影響提供可靠解答。