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

2021-03-05 01:47:24裴志勇吳衛(wèi)國

武漢理工大學學報(交通科學與工程版) 2021年1期

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

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

0 引言

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

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

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

1 結構動力穩(wěn)定性控制方程

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

(1)

式中：U為彈性體系的應變能；V為外力勢能；T為結構的動能.

(2)

可以得到無阻尼情況下結構的平衡方程：

(3)

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

現(xiàn)考慮動載荷為如下形式的周期性載荷，Pcr為結構的靜態(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)，利用正余弦項之間線性獨立的關系，令其系數(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 平面應力四節(jié)點等參元

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

(8)

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

(9)

(10)

則由膜應變對應的單元剛度矩陣Kme為

(11)

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

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

2.2 彎曲單元

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

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

(12)

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

(13)

(14)

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

(15)

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

(16)

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

圖1 MITC等參元和MITC4殼單元

(17)

(18)

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

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

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

(19)

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

扭轉(zhuǎn)虛應變

(20)

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

(21)

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

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

(22)

廣義幾何應變：

(24)

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

(25)

3 計算結果及分析

3.1 加筋板模型

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

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

圖3 加強筋截面尺寸

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

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

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

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

取細長比為γ=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ū).計算結果見圖4，在由β和動載荷頻率θ構成的參數(shù)平面中，隨著靜載分量系數(shù)α的增大，加筋板的動力不穩(wěn)定區(qū)失穩(wěn)頻率值越小，且不穩(wěn)定頻率區(qū)間范圍變得越大，這說明靜載分量越接近臨界載荷，加筋板結構對動載荷的響應越敏感，即使在較低頻率、較小動載荷幅值的動載荷作用下，也可能會發(fā)生失穩(wěn)現(xiàn)象.

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

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

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