結構支撐位置改變時固有頻率的快速計算

歐陽炎，王 棟

(西北工業大學 航空學院 結構動力學與控制研究所，西安 710072)

在結構動力學優化設計過程中，為了獲得滿意的動態響應，需要對結構進行反復修改和設計，其動力學特性的重新分析就不可避免。一般情況下，每次修改后，需先對結構重新進行網格劃分，然后計算結構的固有特性，這是造成結構優化設計時間過長，效率低下的一個主要原因。為了減少有限元重復分析次數，降低優化過程中的計算成本，動響應計算技巧和方法需給予足夠的重視，對于大型復雜結構的動力學問題有必要尋求快速的計算方法［1］。一種行之有效的策略是充分利用優化設計過程中的中間結果，不用直接求解結構修改后的特征方程，根據已得到的結構動態特性，快速、準確地計算模型修改后結構的動態特性，從而有效地減少動力學重分析的次數［2］。

在實際工程設計中，常利用附(增)加支撐(承)的方法來提高結構的剛度和固有頻率，避免產生危險的共振現象。有關梁、薄板結構附加支撐的位置和剛度的優化設計研究，已有許多成果［3－7］。現有的研究成果表明［4－7］:為了提高結構的第一階固有頻率值，附加支撐的最優位置并非總是在有限元網格的節點上，有可能出現在單元的內部。

文獻［7］提出利用插值技術代替網格細化以及特征值求解計算。由于固有頻率對支撐位置的函數畢竟是相當復雜的隱函數，而梁的插值函數僅是位置坐標的3次多項式函數，板的插值函數是位置坐標的(交叉)4次多項式函數，因此人們不禁會對插值計算固有頻率所得結果的精度和適應性產生疑問。本文著重探討支撐位于單元內部時，細長梁和薄板結構固有頻率快速計算方法的精度和適應性問題。

1 支撐剛度等效方法

圖1 單元內部有一個彈性支撐的梁、板單元Fig.1 A beam or plate element with an elastic support

如圖1所示，梁、板單元內部附加有一個彈性支撐。若只考慮支撐的橫向線剛度，并忽略其自身的質量。則支撐的單位等效剛度矩陣Ks只是支撐位置的函數，可由單元的位移形狀函數表示為［6］:

其中，N(a)和N(a，b)分別是梁和板單元的形狀函數矩陣。2節點4自由度梁單元的形狀函數為［9］:

N(ξ)= ［N1N2N3N4］T

其中:

4節點12自由度矩形薄板單元的形狀函數矩陣N(a，b)見文獻［9］。在結構有限元網格不變的情況下，將附加支撐的等效剛度矩陣與結構的剛度矩陣疊加，則系統的頻率特征方程可寫成:

其中，k是支撐的線性剛度系數，Kb、Mb分別為梁或板結構的剛度矩陣和質量矩陣。通過求解特征值問題(3)便可計算出總體結構的固有頻率值ω。可見這樣處理彈性支撐后，求解的特征值問題階數不變。而且Kb和Mb不受支撐位置變化的影響。

“科學是人類認知世界不竭的長河，技術是人類對生存發展方式不倦的創造.研究科學史，本質上也就是研究人類創造的歷史，繼往而開來，有著十分重要的價值和意義.”

2 利用形函數插值方法

對于梁或框架結構，第i階固有頻率ωi相對于支撐位置a的一階導數可表示為［7］:

其中，wi(a)和θi(a)分別為第i階振型在支撐節點處的位移和轉角，而且振型已按質量正交歸一化處理。在計算出彈性支撐位于單元節點上的第i階頻率和對支撐位置的一階導數后，利用梁單元自身的形狀函數N1～N4，可以用插值計算當支撐位于單元內部a處時(如圖1(a)所示)的頻率值:

對于薄板單元，支撐在單元內部的位置如圖1(b)所示。與梁單元類似，可得第i階頻率關于支撐位置的靈敏度為［7］:

可得到支撐位于該單元內部任意點(a，b)處，結構的第i階固有頻率:

其中，wi(a，b)，θxi(a，b)，θyi(a，b)分別為第 i階振型在支撐點(a，b)處的位移、繞x軸轉角和繞y軸轉角。與梁單元類似，在計算出支撐位于單元節點上時的第i階頻率ωi1～ωi4以及頻率關于支撐位置的靈敏度后，利用板單元自身的形函數N1～N4y，并結合形函數在節點上的特性，例如在節點1處有:

利用單元形函數插值法計算固有頻率無需求解特征值問題，僅僅是一些代數運算，因而計算效率更高。由于附加支撐的剛度值和位置會對結構的頻率和振型產生影響，在利用式(5)或式(8)計算固有頻率時，必須考察支撐位于不同節點處時的振型是否一致。

3 數值算列

為了對比兩種方法計算結構固有頻率的精度和效率，將兩種方法用Matlab?編程，分別計算梁和薄板結構附加彈性支撐后的第一階固有頻率值。

3.1 附加彈性支撐的梁結構

如圖2所示懸臂梁，截面為正方形，邊長1 cm，長度L=1 m，密度 ρ=7 800 kg/m3，彈性模量 E=200 GPa。將結構分為10個2節點4自由度梁單元。分別取彈性支撐的剛度系數 k1=3.74 kN/m、k2=16.67 kN/m和k3=44.48 kN/m。根據上述兩種方法，分別計算支撐位置b在梁上連續移動時，懸臂梁第一階固有頻率值，變化曲線如圖3所示。從該圖中可以看出，當彈性支撐在整個梁跨度之間移動時，具有較小剛度值(k1、k2)的支撐采用插值方法計算結果與等效剛度法所得結果基本一致。從計算效率方面來看，支撐正好在網格節點處時進行9次特征值計算后，支撐點位于單元內部時，剛度等效法要再進行多次的特征值求解計算，而插值方法由代數運算即可得出結果。

圖2 附加一個彈性支撐的梁結構Fig.2 A cantilever beam with an elastic support

當支撐的剛度值較大為 k3，并位于區間(0.7，0.8)m時，兩種方法計算結果與理論分析方法［10］所得結果列于表1。從表中可以看出，支撐剛度等效方法與理論解符合較好，計算精度較高。

由文獻［5］中結論結合本文算例的數據可知，當k3=44.48 kN/m，b=0.783 4 m 時，附加彈性支撐可使懸臂梁的第一階頻率值上升到原結構(無支撐時)的第二階頻率。當彈性支撐分別位于相應單元的兩側節點上，即b=0.7 m和b=0.8 m時，懸臂梁的第一、第二階振型如圖4所示，從圖中可以看出，支撐在單元兩側節點處時，第一階振型發生了改變。則根據靈敏度公式(4)計算頻率的一階導數時，對應的振型不同，利用插值公式(5)計算所得結果會有一定的誤差。這是造成圖3中兩種方法所得結果局部不完全相同的原因。由表1可以看出，越靠近單元中間，方法2的誤差越大。

圖3 懸臂梁第一階頻率隨支撐位置b的變化Fig.3 Variations of the 1st frequency of the cantilever beam versus the support position b

圖4 剛度為k3的支撐位于不同節點時懸臂梁的前兩階振型Fig.4 The 1st and 2nd vibration mode shapes of the cantilever beam with an elastic support of stiffness k3located at different nodes

表1 附加彈性支撐懸臂梁的第一階頻率值Tab.1 The 1st frequency of the cantilever beam

3.2 附加彈性支撐的懸臂板結構

如圖5所示正方形懸臂板［7］，邊長 L=W=305 mm，厚度 h=3.28 mm，泊松比 μ =0.3，密度 ρ=2 821 kg/m3，彈性模量E=73.1 GPa。采用4節點12自由度矩形薄板單元，將懸臂板劃分成10×10個單元。

在自由邊上對稱地布置兩個相同剛度的彈性支撐。當支撐數較多時，公式(6)可寫成對多個支撐靈敏度值相加的形式。圖6顯示了不同支撐剛度值時，懸臂板第一階固有頻率隨支撐位置a的變化關系。其中，k1=25.39 kN/m，k2=50.79 kN/m，從圖中可以看出，支撐剛度為k1時，兩種快速計算方法結果基本一致。然而剛度等效方法在處理支撐位于單元內部時，要多次求解高階矩陣(本文為330階)的廣義特征值，因此計算效率很低。

圖5 附加兩個彈性支撐的矩形懸臂板結構Fig.5 A cantilever plate with two elastic supports

當支撐剛度為k2，支撐位置a/W在區間［0.1，0.2］上時，兩種方法所得結果有一定的誤差。當支撐位于 a=30.5 mm(a/W=0.1)，和 a=61 mm(a/W=0.2)時，懸臂板的振型分別如圖7所示。可知懸臂板的第一階振型由扭曲振型轉變為彎曲振型，則由式(6)計算得到的靈敏度所用的振型完全不同，因而由式(8)插值計算得出的頻率值與真實值之間有一定誤差。

在不同的剛度值下，懸臂板的第一階扭曲振型和第一階彎曲振型對應的頻率與支撐位置a的關系如圖8所示。從圖中可以看出，當支撐剛度為k1，a在整個區間上變化時，懸臂板的第一階頻率始終對應第一階彎曲振型，故兩種方法結果吻合很好。當支撐剛度為k2，a在區間［30.5，61］mm時，結構的第一階頻率對應的振型由扭曲振型逐漸變為彎曲振型，即結構的振型在a=30.5 mm和a=61 mm時發生了根本性改變，所以兩種方法的計算結果有誤差，其中支撐剛度等效方法的結果較準確。而當a＞61 mm時，第一階頻率對應振型保持不變，兩種方法計算結果基本一致。

作為對比，若使兩個剛度相同的彈性支撐沿單元的對角線方向移動，如圖5中虛線所示，支撐位置與懸臂板第一階固有頻率關系如圖9所示。其中，k1=25.39 kN/m，k2=76.18 kN/m，從圖中可以看出，支撐剛度為k1時，兩種方法計算結果基本相同。由于插值方法只需在支撐移動軌跡上不同單元的節點處進行少數幾次特征值求解，因此，其計算效率更高。

當支撐剛度為 k2，支撐位置 b/W 在區間［0.1，0.2］，即支撐移動軌跡的第2個單元內時，兩種方法的計算結果有明顯差別。這同樣是由于懸臂板的第一階振型發生改變的結果。其中，支撐剛度等效方法的結果較準確，插值計算結果有一定的誤差。

4 結論

用有限元法計算附有彈性支撐的梁、板結構固有頻率時，如果支撐位置移動到了單元的內部，應用支撐剛度等效或插值方法，都無需重新劃分網格，可以避免因網格變化對計算精度造成的影響。其中，剛度等效法在不增加求解階數的情況下，通過解特征值問題獲得結構的固有頻率。插值法用單元的形函數和固有頻率的導數值計算支撐在單元內部時的頻率值，充分利用了前期的計算結果，因此效率較高。

當彈性支撐的剛度相對較小時，兩種方法的計算結果吻合良好，且都有較高的精度。當支撐剛度較大時，在結構的某些單元節點位置處，結構對應于同一頻率的振型將發生改變，導致頻率關于支撐位置的靈敏度計算公式應用了不同振型的響應值。由此引起利用形函數插值法計算結果會有一定誤差，使用時必須首先加以判斷。支撐剛度等效法無此限制，故應用范圍更廣。

在實際應用中，可充分利用兩種方法各自的優點，將兩種方法結合使用，即在振型可能有突變的單元處，用支撐剛度等效法計算結構的固有頻率，保證計算結果的精度。而在振型無突變的其他單元處，則可考慮用形函數插值方法，以提高計算速度。

［1］ Su H C，Feng R.A new method of structural modal reanalysis for topological modifications［J］.Finite Elements in Analysis and Design，2002，38:1015－1028.

［2］ 何建軍，姜節勝，康興無.結構拓撲大修改的動力學重分析方法［J］.振動工程學報，2007，20(4):407－411.

［3］ 宋 華，胡 瑞，任輝啟，等.內部為彈性點支承的簡支板的振動分析［J］.地下空間與工程學報，2007，3(4):613－616.

［4］ 燕樂維，陳樹輝.一種改進的廣義遺傳算法及其在結構動力優化問題中的應用［J］.工程力學，2010，27(5):21－26.

［5］ Wang C Y.Minimum stiffness of an internal elastic support to maximize the fundamental frequency of a vibrating beam［J］.Journal of Sound and Vibration，2003，259(1):229－232.

［6］ Friswell M I，Wang D.The minimum support stiffness required to raise the fundamental natural frequency of plate structures［J］.Journal of Sound and Vibration，2007，301:665－677.

［7］ Wang D，Jiang J S，Zhang W H.Optimization of support positions to maximize the fundamental frequency of structures［J］. International Journal for Numerical Methods in Engineering，2004，61:1584－1602.

［8］ 李 濤，左正興，廖日東.結構仿真高精度有限元網格劃分方法［J］.機械工程學報，2009，45(6):304－308.

［9］ 朱伯芳.有限單元法原理與應用(第三版)［M］.北京:中國水利水電出版社，2009.

［10］ 方 同，薛 璞.振動理論及應用［M］.西安:西北工業大學出版社，1997.