中厚圓柱殼振型的有限元超收斂拼片恢復(fù)解和網(wǎng)格自適應(yīng)分析

王永亮

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

結(jié)構(gòu)和彈性體的動(dòng)力分析是結(jié)構(gòu)抗震、巖體誘發(fā)地震研究的重要基礎(chǔ)[1-2]，作為支撐結(jié)構(gòu)或儲(chǔ)藏腔體的圓柱殼在結(jié)構(gòu)工程、巖體工程、航空航天工程中廣泛應(yīng)用，研究該結(jié)構(gòu)的振動(dòng)、失穩(wěn)、屈曲等動(dòng)力特性對于研判其失效行為具有重要意義[3-4]。目前研究殼體問題常采用傳統(tǒng)的薄殼理論[5-8]，該理論基于Kirchhoff-Love假定、忽略橫向剪切變形[9]，對于具有較小剪切剛度（即容易發(fā)生顯著橫向剪切變形）的殼體結(jié)構(gòu)將引入一定誤差[10-11]。另外，實(shí)際工程中板殼結(jié)構(gòu)壁厚的增加往往超出薄壁理論的應(yīng)用范圍，也要求考慮橫向剪切變形的影響。中厚殼自由振動(dòng)理論相比薄殼自由振動(dòng)理論，考慮了橫向剪切變形和轉(zhuǎn)動(dòng)慣量的影響，使解答更加可靠。

復(fù)雜結(jié)構(gòu)形式和邊界條件的中厚圓柱殼自由振動(dòng)分析多采用有限元法進(jìn)行求解，高精度自振頻率和振型的數(shù)值解也成為結(jié)構(gòu)精準(zhǔn)損傷識別的關(guān)鍵基礎(chǔ)[12-13]。相比常規(guī)有限元法，自適應(yīng)有限元法可獲得滿足用戶指定誤差限的高精度解答[14]，自適應(yīng)算法成為優(yōu)化網(wǎng)格、提高解答精度的重要方法[15-17]，主要包括提高單元階次的p型自適應(yīng)方法[18]、增加網(wǎng)格密度的h型自適應(yīng)方法[19]、以及二者結(jié)合的ph型自適應(yīng)方法[20]等。自適應(yīng)有限元法在梁、板、殼自由振動(dòng)和彈性穩(wěn)定等結(jié)構(gòu)特征值問題的高精度求解中顯示出良好潛力[21-23]。本文研究中厚圓柱殼的自由振動(dòng)問題，該問題為典型的二階常微分方程組特征值問題；對于各類邊界條件、不同環(huán)向波數(shù)和厚長比中厚圓柱殼的連續(xù)階頻率和振型求解，成為挑戰(zhàn)常微分方程組特征值問題求解的難點(diǎn)。針對上述問題，本文將推廣前期建立的變截面Euler-Bernoulli梁自由振動(dòng)分析的自適應(yīng)算法，合理引入位移超收斂拼片恢復(fù)方法和高階形函數(shù)插值技術(shù)、能量模誤差估計(jì)技術(shù)、網(wǎng)格細(xì)分加密方法，形成一套中厚圓柱殼自由振動(dòng)分析的h型有限元自適應(yīng)求解方案。經(jīng)數(shù)值算例檢驗(yàn)，該方法展示出良好的求解效力，本文旨在報(bào)道這些新進(jìn)展和成果。

1 振動(dòng)方程

考慮圖1所示的中厚圓柱殼，ox為旋轉(zhuǎn)軸，建立中面上A點(diǎn)的局部坐標(biāo)系αβγ，其中α沿子午線切向方向、β沿緯圓切向方向、γ沿法線方向。中厚圓柱殼振動(dòng)的5個(gè)獨(dú)立位移為:沿α方向的線位移u、沿β方向的線位移v、沿γ方向的線位移w、繞α的轉(zhuǎn)角位移?、繞β的轉(zhuǎn)角位移ψ。記圓柱殼中面半徑為r，截面厚度為h，截面剪切剛度修正系數(shù)為κ，慣性矩為J，長度為l；材料彈性模量為E，剪切模量為G，Poisson比為μ，密度為ρ。

圖1 中厚圓柱殼坐標(biāo)系和符號Fig.1 Coordinate systems and symbols of moderately thick circular cylindrical shell

本文研究的中厚圓柱殼自由振動(dòng)的微分控制方程為

式中:（）′=d（）/d x；ω為自振頻率；n為環(huán)向振型的波數(shù)；其余參數(shù)參見文獻(xiàn)[24]。

上述二階常微分方程組特征值問題，可表示為如下矩陣形式

式中:L，R為相應(yīng)的微分算子矩陣；u={u，v，w，?，ψ}T={u1，u2，u3，u4，u5}T為振型（位移）函數(shù)向量；ω，u分別對應(yīng)于特征值和特征向量，本文亦將（ω，u）合稱為特征對。

2 自適應(yīng)目標(biāo)和步驟

設(shè)要求解中厚圓柱殼自由振動(dòng)問題式（2）的某階特征對（ω，u），本文的中厚圓柱殼自由振動(dòng)自適應(yīng)求解目標(biāo)為:在精確解（ω，u）未知的情況下，事先給定誤差限Tol，尋求一個(gè)優(yōu)化的有限元網(wǎng)格π，使該網(wǎng)格上的有限元解（ωh，uh）同時(shí)滿足

實(shí)施時(shí)，由于精確解（ω，u）事先未知，上述目標(biāo)不能作為停機(jī)準(zhǔn)則。對于振型，由于下文介紹的振型（位移）超收斂拼片恢復(fù)解u*較uh具有更高的收斂階，即u*比uh更接近精確解u，因此采用u*代替u對uh進(jìn)行誤差估計(jì)和控制，形成能量模形式的誤差估計(jì)；對于頻率，利用位移超收斂解計(jì)算Rayleigh商[25]得到頻率的超收斂解。通過Rayleigh商得出的頻率超收斂解比振型超收斂解具有更高的收斂階，本文通過估計(jì)和控制振型解答的誤差可保證頻率解答滿足誤差要求，如此形成本方法通過控制振型誤差的停機(jī)準(zhǔn)則?；谟邢拊獾恼`差估計(jì)進(jìn)行網(wǎng)格細(xì)分加密，即可求得優(yōu)化的網(wǎng)格和滿足誤差限的各階高精度頻率和振型解答，形成一套h型有限元自適應(yīng)方案。自適應(yīng)求解步驟包括:

步驟1有限元解——在當(dāng)前網(wǎng)格上（初始網(wǎng)格由用戶提供），進(jìn)行常規(guī)有限元計(jì)算，得到當(dāng)前網(wǎng)格下的有限元解（ωh，uh）；

步驟2誤差估計(jì)——利用有限元后處理超收斂算法，可得當(dāng)前網(wǎng)格下頻率和振型的超收斂解（ω*，u*）；同時(shí)，可得超收斂解與有限元解之間的誤差估計(jì)值；

步驟3網(wǎng)格細(xì)分——對于誤差估計(jì)不滿足誤差限的單元，用網(wǎng)格細(xì)分加密方法將其細(xì)分，得到新的有限元網(wǎng)格，返回步驟1；如果所有單元均滿足誤差限，則網(wǎng)格無需再細(xì)分，求解過程結(jié)束。

3 有限元解

對于求解中厚圓柱殼自由振動(dòng)問題式（2），在給定有限元網(wǎng)格π下，常規(guī)有限元將建立如下線性矩陣特征值問題

式中:D為振型向量；K和M為靜力剛度矩陣和一致質(zhì)量矩陣。在當(dāng)前網(wǎng)格下求解可得有限元解（ωh，uh），需要指出的是，該解答的精度與網(wǎng)格相關(guān)，高精度的解答需要高質(zhì)量的有限元網(wǎng)格。

4 誤差估計(jì)

有限元位移解在單元內(nèi)部具有hm+1階的超收斂性，在單元節(jié)點(diǎn)上具有h2m階的超收斂性，這些節(jié)點(diǎn)成為位移超收斂點(diǎn)。超收斂點(diǎn)上解答的誤差比其他區(qū)域解答更小，具有更快的收斂速度。通過多個(gè)相鄰單元進(jìn)行連接拼片，利用拼片上位于多個(gè)超收斂點(diǎn)的超收斂解和次高階形函數(shù)（高于當(dāng)前形函數(shù)階次m，即≥m+1）插值，即可恢復(fù)獲得拼片單元上的超收斂解，形成位移超收斂拼片恢復(fù)方法[26-27]。本文對于中厚圓柱殼的自由振動(dòng)問題，求得振型（位移）的有限元解后，同樣可以獲得各階振型的超收斂解。通過將單元e的相鄰單元進(jìn)行組合拼片，拼片單元的位移超收斂解的插值形式為

式中:P為給定函數(shù)向量；a為待定系數(shù)向量。

系數(shù)a的取值通過求得如下泛函的極小值使式（6）中乘積結(jié)果在單元節(jié)點(diǎn)處等于當(dāng)前節(jié)點(diǎn)位移值來確定

式中，n為進(jìn)行拼片單元的節(jié)點(diǎn)數(shù)目。

使用最小二乘法求解式（8），即可獲得系數(shù)a的值

式中:各系數(shù)矩陣為

確定a后，即可通過式（6）得到拼片單元的位移超收斂解，該超收斂解具有比當(dāng)前網(wǎng)格下有限元解更高的收斂階，可用于估計(jì)有限元解的誤差。

使用獲得的位移超收斂解，并通過如下Rayleigh商的計(jì)算，可以得到頻率的超收斂解ω*

式中，a（），b（）為應(yīng)變能、動(dòng)能內(nèi)積。利用當(dāng)前網(wǎng)格下位移超收斂解和有限元解，可得能量模形式下的誤差估計(jì)[28]

式中:ne為進(jìn)行拼片的單元數(shù)目；e*為位移超收斂解u*和有限元解uh的誤差，并有能量范數(shù)表達(dá)式

進(jìn)一步，式（12）可寫為

式中，ξ為相對誤差，該值應(yīng)滿足

這里需要指出，通過Rayleigh商計(jì)算得到的頻率超收斂解比振型超收斂解具有更高的收斂階，針對振型進(jìn)行誤差估計(jì)和控制，即可獲得高精度的振型并確保頻率的準(zhǔn)確性。

5 網(wǎng)格細(xì)分加密

如果當(dāng)前單元e上誤差估計(jì)式（15）不滿足，表明該單元需要進(jìn)行細(xì)分加密來降低解答誤差。本文采用將當(dāng)前單元均勻細(xì)分形成多個(gè)單元的網(wǎng)格細(xì)分加密方式，新單元的長度與當(dāng)前單元上的誤差和單元階次相關(guān)，估計(jì)式為

式中:hnew為新單元長度；hold當(dāng)前單元e的長度。為使自適應(yīng)計(jì)算更加高效、避免網(wǎng)格冗余，本文對每次單元細(xì)分加密數(shù)目進(jìn)行控制，實(shí)施方案為

式中:nnew為網(wǎng)格細(xì)分后新單元數(shù)目；d為單元細(xì)分的極限數(shù)目；?·」為向下取整符號。

經(jīng)過以上有限元解的誤差估計(jì)和網(wǎng)格細(xì)分加密，最終可以得到優(yōu)化的網(wǎng)格，并獲得該優(yōu)化網(wǎng)格上滿足誤差限的各階高精度頻率和振型解答。

6 數(shù)值算例

本文算法已經(jīng)編制成Fortran 90程序，本章給出使用本文方法求解多種中厚圓柱殼自由振動(dòng)的數(shù)值算例，通過對比典型的中厚殼、薄殼結(jié)果來驗(yàn)證本文方法計(jì)算結(jié)果的精確性，并討論求解不同環(huán)向波數(shù)、厚長比、邊界條件下中厚圓柱殼自由振動(dòng)問題的適用性。本節(jié)所有算例均采用3次元，初始網(wǎng)格采用2個(gè)單元，給定的初始誤差限為Tol=10-4，設(shè)定單元細(xì)分的極限數(shù)目為d=6。文中k表示階次，置于頻率或振型符號的下標(biāo)位置。

例1薄圓柱殼求解檢驗(yàn)

考慮中厚圓柱殼在兩端簡支條件下的自由振動(dòng)，圓柱殼基本參數(shù)為

采用本文方法求解不同環(huán)向波數(shù)（n=1～15）時(shí)的前2階（k=1，2）特征對，表1所示為計(jì)算得到的頻率結(jié)果。文獻(xiàn)[29]、文獻(xiàn)[30]采用中厚殼理論，分別利用混合有限元法和動(dòng)力剛度法對該問題進(jìn)行了求解，從表1可以看出本文結(jié)果與Sivadas等和陳旭東等的求解結(jié)果展示出良好的吻合性，檢驗(yàn)了本文方法求解各階頻率的可靠性。

表1 中厚圓柱殼自由振動(dòng)頻率值ω/HzTab.1 Frequenciesω/Hz of moderately thick circular cylindrical shell

圖2～圖5列出了n=1和n=10的第1階振型，并在橫坐標(biāo)軸上標(biāo)出自適應(yīng)求解的最終網(wǎng)格分布。需要說明的是，為方便直觀顯示和對比分析，本文中振型結(jié)果均進(jìn)行歸一化處理（令最大振型值為1），同時(shí)在振型圖中也將橫坐標(biāo)軸進(jìn)行歸一化處理（令水平坐標(biāo)軸為x/l）。從圖2、圖3可以看出，對于n=1時(shí)的第1階振型在兩端區(qū)域變化比較劇烈，自適應(yīng)方法劃分出相對比較細(xì)密的網(wǎng)格；振型在中間區(qū)域相對平緩，只需要使用稀疏的網(wǎng)格。

圖2 n=1時(shí)，（，，）振型圖和自適應(yīng)網(wǎng)格Fig.2 Vibration mode（，，）and adaptive mesh with n=1

圖3 n=1時(shí)，（，）振型圖和自適應(yīng)網(wǎng)格Fig.3 Vibration mode（，）and adaptive mesh with n=1

從圖4、圖5可以看出，對于n=10時(shí)的第1階振型在兩全域上變化比較平緩，自適應(yīng)方法劃分出相對稀疏、均勻的網(wǎng)格。有限元網(wǎng)格根據(jù)振型變化情況進(jìn)行自適應(yīng)加密優(yōu)化，展示出良好適用性。

圖4 n=10時(shí)，（，，）振型圖和自適應(yīng)網(wǎng)格Fig.4 Vibration mode（，，）and adaptive mesh with n=10

圖5 n=10時(shí)，（，）振型圖和自適應(yīng)網(wǎng)格Fig.5 Vibration mode（，）and adaptive mesh with n=10

例2中厚圓柱殼連續(xù)階頻率求解

考慮中厚圓柱殼在兩端簡支條件下的自由振動(dòng)，該圓柱殼基本參數(shù)為

使用本文方法計(jì)算了n=0情況下該圓柱殼自由振動(dòng)的前15階（k=1～15）連續(xù)特征對，將結(jié)果列于表2。文獻(xiàn)[31]采用薄殼理論對本圓柱殼進(jìn)行分析，研究中忽略該圓柱殼橫向剪切變形，并使用動(dòng)力剛度法對其進(jìn)行了求解；文獻(xiàn)[32]得到了相應(yīng)的解析形式解答。

表2 中厚圓柱殼自由振動(dòng)頻率值ω/（rad·s-1）Tab.2 Frequenciesω/（rad·s-1）of moderately thick circular cylindrical shell

從表2可以看出本文結(jié)果與陳旭東和Leissa的求解結(jié)果在低階時(shí)吻合很好，展示出本文自適應(yīng)方法對圓柱殼連續(xù)階次特求解的適用性；需要重點(diǎn)說明的是，隨著階次升高（如k=11～15），頻率結(jié)果的差異愈加明顯，這是由于采用薄殼理論對中厚圓柱殼結(jié)構(gòu)進(jìn)行分析引入的誤差。

例3中厚圓柱殼不同環(huán)向波數(shù)n和厚長比h/l

下面分別考慮兩端簡支中厚圓柱殼在不同環(huán)向波數(shù)n、厚長比h/l條件下自由振動(dòng)的求解。首先，考慮中厚圓柱殼在不同環(huán)向波數(shù)n（n=1～5）時(shí)的自由振動(dòng)，設(shè)定圓柱殼Poisson比μ=0.3、長徑比l/r=1.0、厚徑比h/r=0.1，其余參數(shù)為

采用本文方法求解各種工況下的第1階特征對，將自振頻率結(jié)果轉(zhuǎn)化為無量綱值列于表3，文獻(xiàn)[33]采用沿圓柱殼周向位移分布形狀函數(shù)的有限元模型、文獻(xiàn)[34]采用層狀圓柱殼模型對本算例求解并得到有效的解答，從表3可以看出的本文方法求得不同環(huán)向波數(shù)下中厚圓柱殼圓柱殼的結(jié)果與文獻(xiàn)解答吻合性較好，展示出本文自適應(yīng)方法對于各類環(huán)向波數(shù)圓柱殼體的適用性。

表3 中厚圓柱殼不同環(huán)向波數(shù)n自由振動(dòng)無量綱頻率值Tab.3 Non-dimensional frequencies of moderately thick circular cylindrical shell with different circumferential wave number n

表3 中厚圓柱殼不同環(huán)向波數(shù)n自由振動(dòng)無量綱頻率值Tab.3 Non-dimensional frequencies of moderately thick circular cylindrical shell with different circumferential wave number n

環(huán)向波數(shù)n頻率值ω-h1 ω-h1[33] ω-h1[34]1 1.063 43 1.062 26 1.062 34 2 0.881 21 0.882 33 0.882 53 3 0.806 67 0.809 25 0.809 51 4 0.897 56 0.898 77 0.898 93 5 1.115 64 —1.122 09

同時(shí)，考慮中厚圓柱殼在不同厚長比h/l（h/l=0.01，h/l=0.10，h/l=0.20，h/l=0.40，h/l=0.60，h/l=0.80，h/l=1.00）工況下自由振動(dòng)，設(shè)定圓柱殼環(huán)向波數(shù)n=0、徑厚比h/r=0.4，采用本文方法求得了不同厚長h/l比下的第1階（k=1）特征對，并將自振頻率結(jié)果轉(zhuǎn)化為無量綱值列于表4，可以看出本文結(jié)果與Armenakas等和Loy等求解的近似結(jié)果吻合較好，展示出本文自適應(yīng)方法對于薄壁、厚壁等不同厚長比圓柱殼體的適用性。

表4 中厚圓柱殼不同厚長比h/l自由振動(dòng)無量綱頻率值Tab.4 Non-dimensional frequencies of moderately thick circular cylindrical shell with different ratio of thickness to length h/l

表4 中厚圓柱殼不同厚長比h/l自由振動(dòng)無量綱頻率值Tab.4 Non-dimensional frequencies of moderately thick circular cylindrical shell with different ratio of thickness to length h/l

厚長比h/l頻率值ω-h1 ω-h1[33] ω-h1[34]0.01 0.016 11 0.016 12 0.010 02 0.10 0.152 91 0.152 89 0.100 00 0.20 0.203 76 0.204 95 0.200 00 0.40 0.278 54 0.275 40 0.275 44 0.60 0.422 12 0.420 22 0.420 35 0.80 0.597 75 0.600 09 0.600 33 1.00 0.784 32 0.792 74 0.793 14

例4中厚圓柱殼不同邊界條件

為檢驗(yàn)本文方法對不同邊界條件分析的適用性，下面分別對兩端簡支、兩端固定等典型邊界條件的中厚圓柱殼自由振動(dòng)進(jìn)行求解。設(shè)中厚圓柱殼Poisson比μ=0.3、厚徑比h/r=0.4，其余參數(shù)為

采用本文方法求解在兩端簡支、兩端固定邊界時(shí)圓柱殼環(huán)向波數(shù)n=0、不同厚長比h/l（h/l=0.1，h/l=0.2，h/l=0.4，h/l=0.6，h/l=0.8，h/l=1.0）工況下第2階（k=2）特征對，將自振頻率結(jié)果轉(zhuǎn)化為無量綱值列于表5，可以看出本文結(jié)果與Loy等采用層狀圓柱殼模型求解結(jié)果相吻合，展示出本文自適應(yīng)方法對不同邊界條件圓柱殼體分析的適用性。需要重點(diǎn)指出的是，本文自適應(yīng)方法兼具通用性，誤差估計(jì)和網(wǎng)格劃分技術(shù)不依賴于邊界條件，有潛力推廣到工程實(shí)際問題中力邊界、彈性邊界、復(fù)合邊界等復(fù)雜邊界工況。

表5 中厚圓柱殼不同邊界條件下自由振動(dòng)無量綱頻率值Tab.5 Non-dimensional frequencies of moderately thick circular cylindrical shell under different boundary conditions

表5 中厚圓柱殼不同邊界條件下自由振動(dòng)無量綱頻率值Tab.5 Non-dimensional frequencies of moderately thick circular cylindrical shell under different boundary conditions

厚長比h/l頻率值兩端簡支ω-h2 ω-h2[34]兩端固支ω-h2 ω-h2[34]0.1 0.932 87 0.929 30 1.026 52 1.035 67 0.2 1.076 53 1.069 48 1.437 37 1.477 66 0.4 2.051 23 2.047 18 3.152 65 3.277 61 0.6 3.595 71 3.621 23 5.099 47 5.323 09 0.8 5.327 11 5.398 23 7.464 93 7.391 31 1.0 7.116 82 7.243 17 9.326 62 9.461 18

圖6 兩端簡支邊界條件下（）振型圖和自適應(yīng)網(wǎng)格Fig.6 Vibration mode（）and adaptive mesh under simply supported boundary conditions at both ends

圖7 兩端固支邊界條件下（，，）振型圖和自適應(yīng)網(wǎng)格Fig.7 Vibration mode（，，）and adaptive mesh under fixed boundary conditions at both ends

7 結(jié) 論

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

（1）本文建立了中厚圓柱殼振型的超收斂拼片恢復(fù)解，利用超收斂解估計(jì)當(dāng)前網(wǎng)格下有限元解的誤差，指導(dǎo)單元加密細(xì)分，形成非均勻分布的優(yōu)化網(wǎng)格。在振型解答變化相對比較劇烈的區(qū)域，自適應(yīng)方法可自動(dòng)劃分出相對細(xì)密的有限元網(wǎng)格，其他區(qū)域則使用相對稀疏網(wǎng)格。

（2）本文方法在中厚圓柱殼簡支、固支等多類邊界條件、不同環(huán)向波數(shù)和厚長比工況下連續(xù)階的頻率和振型求解中，均成功獲得了優(yōu)化的有限元網(wǎng)格和滿足預(yù)設(shè)誤差限的高精度解答，展示出良好的求解效力。

（3）本文方法兼具通用性，有潛力推廣到一般結(jié)構(gòu)特值問題振型（位移場）和固體應(yīng)力（位移導(dǎo)數(shù)場）的精細(xì)化數(shù)值模型和高精度計(jì)算領(lǐng)域，成為下一階段的研究內(nèi)容。