基于MITC單元的加筋板動力穩(wěn)定性分析*

李 念 裴志勇 吳衛(wèi)國

(武漢理工大學(xué)交通學(xué)院1) 武漢 430063)(武漢理工大學(xué)綠色智能江海直達(dá)船舶與郵輪游艇研究中心2) 武漢 430063)

0 引 言

為了確保結(jié)構(gòu)具有足夠的穩(wěn)定性，臨界載荷的計算非常重要.當(dāng)所考慮的結(jié)構(gòu)承受的載荷是動載荷時，則屬于動力穩(wěn)定性問題，在該問題上，由于載荷隨著時間而變化，使得動力穩(wěn)定性問題要比靜力的情況復(fù)雜很多.

在彈性范圍內(nèi)，通常將動力穩(wěn)定性問題轉(zhuǎn)化為一個Mathiue方程進(jìn)行求解，有限元方法也應(yīng)用到動力穩(wěn)定性問題中，Thana等[1]采用有限元法分析討論了梁結(jié)構(gòu)在軸向周期荷載作用下的動力穩(wěn)定性問題.Loja等[2]研究了可變剛度的復(fù)合材料板的動穩(wěn)性，并考慮了靜態(tài)和動態(tài)載荷因數(shù)，以及邊界條件的影響.Talimian等[3]采用有限差分法求解板的運(yùn)動微分方程，分別從載荷系數(shù)、載荷類型，及其組合形式等方面進(jìn)了討論，總結(jié)了載荷參數(shù)對不穩(wěn)定區(qū)的影響規(guī)律.相比于平板而言，加筋板由于其板和筋之間相互耦合，使得從解析角度難以求解.目前的研究方法，一類是將其簡化為便于處理的計算模型，如正交異性板模型、交叉梁系模型[4-5]，但它們在使用上都具有一定的局限性；另一類方法則是對加筋板進(jìn)行離散化處理，采用有限元法、有限差分法等數(shù)值方法進(jìn)行求解，劉媛[6]基于Timoshenko梁和Mindlin板理論，建立了加筋板動力穩(wěn)定性分析的有限元模型，具體分析了加強(qiáng)筋尺寸、種類及其數(shù)目對結(jié)構(gòu)動力穩(wěn)定性的影響.

MITC單元是一種基于假定應(yīng)變法的張量分量混合插值板單元(mixed interpolation of tensorial components，MITC)，它具有完備的理論基礎(chǔ)，不會出現(xiàn)剪切鎖閉現(xiàn)象，且對單元的幾何扭曲不敏感，能很好地適用于各類板殼結(jié)構(gòu)，在收斂性和通用性方面具備優(yōu)勢[7-8].為了探究在工程上受到廣泛應(yīng)用的加筋板結(jié)構(gòu)的動力穩(wěn)定性，文中根據(jù)拉格朗日方程及哈密爾頓原理得到結(jié)構(gòu)動穩(wěn)性的矩陣方程，并利用Bolotin方法將動力穩(wěn)定性問題轉(zhuǎn)換為求解廣義特征值問題，采用MITC4殼單元對加筋板進(jìn)行離散化，推導(dǎo)了四節(jié)點(diǎn)24自由度殼單元特性矩陣，并在Matlab平臺上開發(fā)了求解加筋板屈曲、模態(tài)以及動力穩(wěn)定性的有限元計算程序，以周期性邊界條件下加筋板為研究對象進(jìn)行計算分析，先通過對比MITC4單元以及Abaqus軟件的屈曲和模態(tài)計算結(jié)果，驗證了MITC4殼單元的通用性和準(zhǔn)確性，并能較好地預(yù)報加筋板的動力穩(wěn)定性，最后考慮了靜載系數(shù)、截面參數(shù)等因素對加筋板的動力穩(wěn)定性的影響，為結(jié)構(gòu)穩(wěn)定性的改善提供設(shè)計參考.

1 結(jié)構(gòu)動力穩(wěn)定性控制方程

根據(jù)拉格朗日方程L=U+V-T及哈密爾頓原理，有如下關(guān)系.

(1)

式中：U為彈性體系的應(yīng)變能；V為外力勢能；T為結(jié)構(gòu)的動能.

(2)

可以得到無阻尼情況下結(jié)構(gòu)的平衡方程：

(3)

式中：M為質(zhì)量矩陣；K為剛度矩陣；Kg為幾何剛度矩陣.

現(xiàn)考慮動載荷為如下形式的周期性載荷，Pcr為結(jié)構(gòu)的靜態(tài)臨界載荷.

P(t)=Ps+Pdcos (θt)=Pcr(α+βcos (θt))

(4)

根據(jù)Bolotin方法，形如式(3)帶有周期系數(shù)的二階Mathieu微分方程解的收斂性可以被周期為T和2T(T=2π/θ)的解為邊界所確定，分別將兩類解以傅里葉級數(shù)的形式展開.

(5)

將式(4)～(5)帶入式(3)，利用正余弦項之間線性獨(dú)立的關(guān)系，令其系數(shù)為零可以得到以下線性方程組.

(6b)

(6c)

(6d)

方程組(6)存在非零解的充要條件是其系數(shù)矩陣為奇異矩陣，通過求解前n階行列式可以得到前n個動力不穩(wěn)定區(qū)域，隨著階數(shù)的增大，頻率區(qū)間迅速減小，其中第一階不穩(wěn)定區(qū)已包括參數(shù)平面的絕大部分，也稱為主要動力不穩(wěn)定區(qū)，它可以由式(6a)與(6b)的一階解確定.

(7)

通過求解式(7)矩陣的廣義特征值，從而可以得到由α，β，θ等參數(shù)決定的動力不穩(wěn)定區(qū)域.

2 MITC4殼單元特性矩陣

2.1 平面應(yīng)力四節(jié)點(diǎn)等參元

單元局部坐標(biāo)系下，雙線性四邊形四節(jié)點(diǎn)等參元插值函數(shù)為

(8)

根據(jù)平面應(yīng)力問題的應(yīng)變-位移幾何關(guān)系式，有：

(9)

(10)

則由膜應(yīng)變對應(yīng)的單元剛度矩陣Kme為

(11)

式中：ξm,ηk為等參元在自然坐標(biāo)系下高斯積分點(diǎn)坐標(biāo)；Wm,Wk為相應(yīng)的高斯積分點(diǎn)權(quán)重系數(shù)；J為雅可比矩陣.

數(shù)值積分方法選用高斯積分，雙線性四邊形單元采用2×2階高斯積分即可達(dá)到完全積分效果.

2.2 彎曲單元

目前常用的板彎曲殼單元是基于Mindlin板理論的板單元，這種考慮了板橫向剪切應(yīng)變的單元適合模擬厚板彎曲.但在有限元模擬時，如果劃分單元的厚度與邊長之比較小時，會導(dǎo)致單元的橫向剪切應(yīng)變被夸大，剪切剛度人為“剛化”，即發(fā)生剪切鎖閉現(xiàn)象.盡管采用減縮積分或者選擇積分等方式能解決剪切鎖閉現(xiàn)象，但同時有可能會導(dǎo)致多余的零能模式，使得求解失真.因而，選用一種通用性更好，計算更為精確的單元是有必要的，MITC殼單元能較好地解決剪切自鎖現(xiàn)象，且不必對數(shù)值積分方案做調(diào)整.

彎曲單元位移插值函數(shù)采用與3.1節(jié)平面應(yīng)力單元相同的雙線性函數(shù)：

(12)

由幾何方程得到板面內(nèi)彎曲應(yīng)變關(guān)系：

(13)

(14)

由高斯積分得到單元彎曲剛度矩陣Kbe：

(15)

考慮板橫向剪切應(yīng)變的影響，基于Mindlin板理論，板單元的剪切應(yīng)變?yōu)槭?16)，結(jié)合插值函數(shù)式(12)可得到應(yīng)變矩陣Bs.

(16)

而基于MITC板理論，為了消除剪切閉鎖，使用與Mindlin板單元相同的插值函數(shù)，用插值點(diǎn)處的剪切應(yīng)變值，對插值點(diǎn)處的橫向位移、轉(zhuǎn)角等協(xié)調(diào)應(yīng)變進(jìn)行混合插值，從而用假定的剪切應(yīng)變來替換單元剪切應(yīng)變部分，達(dá)到消除剪切自鎖的目的.以單元四邊中點(diǎn)的剪切應(yīng)變進(jìn)行插值，得到一組假定應(yīng)變，見圖1.

圖1 MITC等參元和MITC4殼單元

(17)

(18)

2.3 第六自由度扭轉(zhuǎn)單元

根據(jù)殼單元的公式，每個節(jié)點(diǎn)的剛度矩陣只耦合5個自由度，然而，采用5自由度殼單元離散加筋板結(jié)構(gòu)存在困難.為此，文中基于Kanok[9]的方法，引入了面內(nèi)扭轉(zhuǎn)自由度分量.

定義罰能量函數(shù)，當(dāng)殼體發(fā)生剛體轉(zhuǎn)動，罰能量變?yōu)榱悖瑸?/p>

(19)

式中：G為材料剪切模量；κT為罰因子，文獻(xiàn)[9]推薦取大于0.1的常數(shù)來抑制面內(nèi)扭轉(zhuǎn)剛度，隨著κT增大，面內(nèi)扭轉(zhuǎn)模態(tài)的影響會迅速減小并趨于穩(wěn)定，文中取κT=1.

扭轉(zhuǎn)虛應(yīng)變

(20)

采用減縮高斯積分計算單元扭轉(zhuǎn)剛度矩陣Kte

(21)

2.4 單元質(zhì)量矩陣與幾何剛度矩陣

單元質(zhì)量矩陣的推導(dǎo)類似于剛度矩陣，通過Hamilton原理即可得到單元質(zhì)量矩陣的一般表達(dá)式.通常，質(zhì)量矩陣分為一致質(zhì)量陣和集中質(zhì)量陣，考慮到動力穩(wěn)定性問題會涉及到頻率計算，文中采用在頻率計算方面更為精確的集中質(zhì)量矩陣.

(22)

廣義幾何應(yīng)變：

(24)

由虛功原理得到幾何剛度矩陣：

(25)

3 計算結(jié)果及分析

3.1 加筋板模型

現(xiàn)考慮三跨加筋板模型，圖2.a=2 550 mm，b=850 mm，材料的彈性模量E=2.06×1011Pa，泊松比ν=0.3，密度ρ=7 850 kg/m3，受到縱向周期載荷的作用.約束條件為周期性邊界條件，即在加筋板的四條邊中，對邊分別具有相同的位移模式，加載端和其中一條鄰邊用MPC約束，加筋板是構(gòu)成大型結(jié)構(gòu)的基本單元，采用周期性邊界條件能較好地反映周邊的影響，從而能夠用部分的性質(zhì)來表達(dá)整體的性質(zhì).加強(qiáng)筋尺寸見圖3，h×bf×tw/tf=600 mm×150 mm×15/20 mm，板和加強(qiáng)筋均用殼單元進(jìn)行離散，網(wǎng)格大小為50 mm.

圖2 加筋板模型的尺寸和邊界條件

圖3 加強(qiáng)筋截面尺寸

3.2 屈曲及模態(tài)分析

表1 屈曲及模態(tài)(一階)對比

3.3 動力穩(wěn)定性分析

3.3.1靜載系數(shù)的影響

取細(xì)長比為γ=2.21的加筋板分析計算，由3.2節(jié)的屈曲分析得到該加筋板的靜態(tài)臨界力Pcr=4 357 N/mm，動載荷P(t)=Ps+Pdcos (θt)=Pcr(α+βcosθt)，其中靜載荷分量Ps=αPcr，為探討靜載分量對加筋板動力穩(wěn)定性的影響，令α=0，0.2，0.4，0.6，計算不同靜載荷分量下加筋板的動力不穩(wěn)定區(qū).計算結(jié)果見圖4，在由β和動載荷頻率θ構(gòu)成的參數(shù)平面中，隨著靜載分量系數(shù)α的增大，加筋板的動力不穩(wěn)定區(qū)失穩(wěn)頻率值越小，且不穩(wěn)定頻率區(qū)間范圍變得越大，這說明靜載分量越接近臨界載荷，加筋板結(jié)構(gòu)對動載荷的響應(yīng)越敏感，即使在較低頻率、較小動載荷幅值的動載荷作用下，也可能會發(fā)生失穩(wěn)現(xiàn)象.

圖4 不同靜載分量下的動力不穩(wěn)定區(qū)域

3.3.2加筋板截面參數(shù)的影響

最小剖面模數(shù)是衡量加筋板強(qiáng)度的重要指標(biāo)，為研究加筋板的截面參數(shù)對動穩(wěn)性的影響，設(shè)計了6組最小剖面模數(shù)基本相同，而細(xì)長比各不相同的加筋板，見表2.分別計算這6組加筋板的動力不穩(wěn)定區(qū)(周期載荷靜載分量均取α=0)，計算結(jié)果見圖5，細(xì)長比是反應(yīng)加筋板柔度的物理量，隨著加筋板柔度的增大，失穩(wěn)頻率逐漸降低，同時不穩(wěn)定區(qū)域的范圍也逐漸縮小.因此，在最小剖面模數(shù)相同的前提下，加筋板的設(shè)計需要綜合考慮外載荷的幅值和頻率范圍，選取合適的尺寸，從而提高加筋板的動力穩(wěn)定性.

表2 加筋板截面參數(shù)

圖5 不同細(xì)長比的動力不穩(wěn)定區(qū)域

4 結(jié) 論

1) MITC4單元在屈曲分析和模態(tài)分析中表現(xiàn)良好，對不同細(xì)長比的加筋板的計算結(jié)果均較好，能很好的模擬薄殼和厚殼，具備良好的通用性，適用于動力穩(wěn)定性計算.

2) 周期載荷的靜載分量對動力穩(wěn)定性有顯著的影響，當(dāng)靜載分量較大時，結(jié)構(gòu)失穩(wěn)頻率降低，不穩(wěn)定區(qū)域范圍也迅速增大，即使在外載荷幅值較小的情況下，也可能發(fā)生失穩(wěn)現(xiàn)象.

3) 在加筋板截面的最小剖面模數(shù)一定的前提下，細(xì)長比越大，失穩(wěn)頻率越小，不穩(wěn)定區(qū)的范圍也減小.加筋板結(jié)構(gòu)設(shè)計時需綜合考慮周期載荷的幅值和頻率，選取合適的截面尺寸以改善動力穩(wěn)定性.