中厚圓柱殼振型的有限元超收斂拼片恢復解和網格自適應分析

王永亮

（1.中國礦業大學（北京） 力學與建筑工程學院，北京 100083；2.中國礦業大學（北京） 煤炭資源與安全開采國家重點實驗室，北京 100083）

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

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

1 振動方程

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

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

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

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

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

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

2 自適應目標和步驟

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

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

步驟1有限元解——在當前網格上（初始網格由用戶提供），進行常規有限元計算，得到當前網格下的有限元解（ωh，uh）；

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

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

3 有限元解

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

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

4 誤差估計

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

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

系數a的取值通過求得如下泛函的極小值使式（6）中乘積結果在單元節點處等于當前節點位移值來確定

式中，n為進行拼片單元的節點數目。

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

式中:各系數矩陣為

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

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

式中，a（），b（）為應變能、動能內積。利用當前網格下位移超收斂解和有限元解，可得能量模形式下的誤差估計[28]

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

進一步，式（12）可寫為

式中，ξ為相對誤差，該值應滿足

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

5 網格細分加密

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

式中:hnew為新單元長度；hold當前單元e的長度。為使自適應計算更加高效、避免網格冗余，本文對每次單元細分加密數目進行控制，實施方案為

式中:nnew為網格細分后新單元數目；d為單元細分的極限數目；?·」為向下取整符號。

經過以上有限元解的誤差估計和網格細分加密，最終可以得到優化的網格，并獲得該優化網格上滿足誤差限的各階高精度頻率和振型解答。

6 數值算例

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

例1薄圓柱殼求解檢驗

考慮中厚圓柱殼在兩端簡支條件下的自由振動，圓柱殼基本參數為

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

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

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

圖2 n=1時，（，，）振型圖和自適應網格Fig.2 Vibration mode（，，）and adaptive mesh with n=1

圖3 n=1時，（，）振型圖和自適應網格Fig.3 Vibration mode（，）and adaptive mesh with n=1

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

圖4 n=10時，（，，）振型圖和自適應網格Fig.4 Vibration mode（，，）and adaptive mesh with n=10

圖5 n=10時，（，）振型圖和自適應網格Fig.5 Vibration mode（，）and adaptive mesh with n=10

例2中厚圓柱殼連續階頻率求解

考慮中厚圓柱殼在兩端簡支條件下的自由振動，該圓柱殼基本參數為

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

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

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

例3中厚圓柱殼不同環向波數n和厚長比h/l

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

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

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

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

環向波數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

同時，考慮中厚圓柱殼在不同厚長比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）工況下自由振動，設定圓柱殼環向波數n=0、徑厚比h/r=0.4，采用本文方法求得了不同厚長h/l比下的第1階（k=1）特征對，并將自振頻率結果轉化為無量綱值列于表4，可以看出本文結果與Armenakas等和Loy等求解的近似結果吻合較好，展示出本文自適應方法對于薄壁、厚壁等不同厚長比圓柱殼體的適用性。

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

表4 中厚圓柱殼不同厚長比h/l自由振動無量綱頻率值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中厚圓柱殼不同邊界條件

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

采用本文方法求解在兩端簡支、兩端固定邊界時圓柱殼環向波數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）特征對，將自振頻率結果轉化為無量綱值列于表5，可以看出本文結果與Loy等采用層狀圓柱殼模型求解結果相吻合，展示出本文自適應方法對不同邊界條件圓柱殼體分析的適用性。需要重點指出的是，本文自適應方法兼具通用性，誤差估計和網格劃分技術不依賴于邊界條件，有潛力推廣到工程實際問題中力邊界、彈性邊界、復合邊界等復雜邊界工況。

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

表5 中厚圓柱殼不同邊界條件下自由振動無量綱頻率值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 兩端簡支邊界條件下（）振型圖和自適應網格Fig.6 Vibration mode（）and adaptive mesh under simply supported boundary conditions at both ends

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

7 結 論

本文研究工作的主要結論和計算方法潛力可總結如下:

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

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

（3）本文方法兼具通用性，有潛力推廣到一般結構特值問題振型（位移場）和固體應力（位移導數場）的精細化數值模型和高精度計算領域，成為下一階段的研究內容。