王驍峰，段 毅，袁銳之

(空間物理重點(diǎn)試驗(yàn)室， 北京 100076)

硬殼結(jié)構(gòu)是最簡單的薄壁結(jié)構(gòu)[1]，圓柱筒殼由于特殊的幾何構(gòu)型和優(yōu)良的受力性能，被廣泛應(yīng)用于船舶工程、導(dǎo)彈火箭外殼及航空航天等許多工程技術(shù)領(lǐng)域[2]，如火箭、導(dǎo)彈的儀器艙、貯箱、發(fā)動機(jī)燃燒室等。圓筒硬殼結(jié)構(gòu)在火箭等飛行器的飛行過程中，飛行器艙段的薄壁圓筒殼受到軸向壓力往往會在強(qiáng)度破壞之前失穩(wěn)(發(fā)生屈曲)，使飛行器的結(jié)構(gòu)承載失效。均勻軸壓下圓柱薄殼的極限承載能力對初始幾何缺陷非常敏感。Tennyson[3]等人研究了局部軸對稱缺陷對圓柱薄殼屈曲性能的影響，指出局部缺陷會大大降低結(jié)構(gòu)的承載性能。葉軍等的研究也表明[4]，環(huán)向初始缺陷是影響薄壁鋼桶臨界屈曲載荷的主要因素。

對于圓筒硬殼結(jié)構(gòu)來說，由于加工工藝、加工誤差等原因，使得圓筒殼不可避免的存在幾何缺陷，即存在縱向和環(huán)向結(jié)構(gòu)上的不對稱性，如焊縫、同軸度誤差、不圓度誤差以及局部凹凸不平等等，這種不對稱性導(dǎo)致臨界軸壓明顯降低。造成缺陷的原因還有存儲、運(yùn)輸、安裝和使用等[5]。但加工誤差所造成的初始幾何缺陷是殼體承載能力下降的主要原因[6]。在進(jìn)行殼體屈曲特性研究時，最為理想的情況是引入真實(shí)幾何缺陷，但是，由于殼體加工工藝不同，批次不同，產(chǎn)生的幾何缺陷也不同，即使同一批次產(chǎn)品的幾何缺陷也具有隨機(jī)性，這就需要通過大批量、相似殼體試驗(yàn)，建立真實(shí)幾何缺陷的可靠性模型，往往費(fèi)時費(fèi)力，收效甚微[6]。因此，引入等效幾何缺陷的方法研究殼體屈曲特性就十分必要。Tennyson[3]研究局部軸對稱缺陷對鋼筒屈曲性能的影響時，發(fā)現(xiàn)即使只有一條軸對稱缺陷，也會使薄壁筒殼的屈曲臨界載荷有較大的降低，Koiter[7]首次研究了非完善殼體穩(wěn)定性的一般準(zhǔn)則，提出了“缺陷敏感度”的概念，揭示了屈曲的跳躍性與原始初始缺陷理論的相互關(guān)系。文獻(xiàn)[8]表明圓筒環(huán)向焊縫缺陷會導(dǎo)致最低的屈曲承載力，而缺陷最不利位置可通過無缺陷圓筒非線性屈曲模態(tài)確定[8]。歐洲法規(guī)EN1993-1-6(2007)[9]規(guī)定，幾何缺陷的形狀應(yīng)取最差缺陷(即導(dǎo)致殼體屈曲載荷下降最快的缺陷)，在最差缺陷形狀未知的情況下，建議采用模態(tài)缺陷分析殼體的屈曲特性。

初始缺陷的形式一般分為特征值屈曲模態(tài)缺陷、環(huán)向焊縫缺陷、隨機(jī)缺陷、實(shí)測缺陷等[10]。對于結(jié)構(gòu)的實(shí)際缺陷來說，在產(chǎn)品失效之前，無法預(yù)知缺陷的性質(zhì)、位置和和缺陷程度，無法對后三種缺陷的屈曲問題進(jìn)行前瞻性分析。工程上常常采用特征值屈曲模態(tài)缺陷進(jìn)行分析，無需知道缺陷的性質(zhì)、位置和和缺陷程度。對于殼體來說，只需給出殼體的線性特征值屈曲模態(tài)和一定百分比的殼體厚度作為初始幾何缺陷。對于網(wǎng)格加筋的圓筒結(jié)構(gòu)、錐形圓筒結(jié)構(gòu)等，求解殼體穩(wěn)定性問題時，均可以當(dāng)量成光壁圓筒結(jié)構(gòu)[11]。本文研究表明，通過基于屈曲模態(tài)缺陷的非線性屈曲分析來求解圓筒結(jié)構(gòu)的臨界失穩(wěn)軸壓，用于工程設(shè)計是可行的。

1 非線性屈曲的描述

對于板、殼、桿、梁等結(jié)構(gòu)，當(dāng)其所受的外載荷達(dá)到一定值的時候，結(jié)構(gòu)失穩(wěn)，此時無需增大外載荷，甚至在外載荷減小的情況下，結(jié)構(gòu)的變形繼續(xù)增大，即結(jié)構(gòu)的剛度為零，甚至為負(fù)。隨著變形的增大，結(jié)構(gòu)又開始恢復(fù)抵抗變形的能力，即結(jié)構(gòu)發(fā)生了屈曲。第一次失穩(wěn)發(fā)生前的過程稱為前屈曲，第一次失穩(wěn)的后繼過程稱為后屈曲。

當(dāng)結(jié)構(gòu)處于理想狀態(tài)、無任何物理缺陷和幾何缺陷時，結(jié)構(gòu)的彈性失穩(wěn)表現(xiàn)為歧點(diǎn)(分岔點(diǎn))失穩(wěn)；對結(jié)構(gòu)的屈曲分析可以用線性特征值描述，通過求解特征值可以確定結(jié)構(gòu)的臨界失穩(wěn)載荷。這個分析過程也是一個典型的前屈曲分析過程，如圖1中理想加載路徑所示。但實(shí)際結(jié)構(gòu)不可能是理想結(jié)構(gòu)，幾何不對稱因素(幾何缺陷)是不可避免的。結(jié)構(gòu)的穩(wěn)定性對初始幾何缺陷十分敏感，結(jié)構(gòu)的失穩(wěn)過程如圖1中所示非理想結(jié)構(gòu)加載路徑。結(jié)構(gòu)失穩(wěn)的形式是駐點(diǎn)失穩(wěn)，失穩(wěn)的過程是一個幾何非線性過程，此時對屈曲過程的描述應(yīng)采用非線性描述，由圖1可以看出，特征值屈曲的歧點(diǎn)所對應(yīng)的臨界壓力值是非保守解，而結(jié)構(gòu)的實(shí)際臨界壓力值應(yīng)該是駐點(diǎn)所對應(yīng)的臨界壓力值。

圖1 結(jié)構(gòu)屈曲的載荷(F)-位移(u)示意圖

2 特征值屈曲

結(jié)構(gòu)在外載荷作用下任意狀態(tài)時的增量平衡方程是

(1)

(2)

達(dá)到失穩(wěn)狀態(tài)時，結(jié)構(gòu)失去進(jìn)一步抵抗變形的能力，此時應(yīng)有ΔPN≈0，而增量平衡方程變?yōu)?/p>

(3)

因此，線性屈曲轉(zhuǎn)化為特征值問題。最小特征值λ代表了臨界載荷比例因子，特征向量ΔuM代表了失穩(wěn)模態(tài)。

3 廣義弧長法求解非線性屈曲

非線性屈曲的求解一般采用廣義弧長法，即每一個增量步是弧長步而不是載荷步。在圖2所示的求解空間里，弧長法是借助一個弧面將載荷因子增量Δλ和位移增量ΔuN關(guān)聯(lián)起來。

其中，弧長由下式?jīng)Q定

(4)

圖2表示一個弧長步內(nèi)的迭代求解過程；一個非線性問題的求解過程是由多個弧長步組成。

4 修正的Riks算法

4.1 初始載荷因子增量的求解

修正的Riks算法即修正的弧長法，其特點(diǎn)是用“平面”代替“弧面”，具有求解精度好、效率高等特點(diǎn)，為ABAQUS所采用。在多維求解空間里，修正的Riks法(Modified Riks algorithm)的求解過程如圖3所示。

圖2 弧長法求解非線性屈曲

圖3 修正的Riks算法求解非線性屈曲

(5)

第一個弧長Δl的大小由用戶指定。則應(yīng)有

(6)

(7)

從而

(8)

此處

(9)

求得Δλ0以后，開始進(jìn)行如下過程的迭代。

載荷因子增量及位移增量初始化，有

Δλi=Δλ0

(10)

(11)

4.2 迭代過程

在第i個迭代步，形成節(jié)點(diǎn)應(yīng)力矩陣IN，切線剛度矩陣KNM

(12)

(13)

其中βN是應(yīng)變位移轉(zhuǎn)換矩陣。這樣，求解空間中的力和位移的預(yù)測點(diǎn)Ai所對應(yīng)的狀態(tài)為

(14)

檢查迭代平衡

(15)

式(15)為解空間點(diǎn)Ai所對應(yīng)的力和節(jié)點(diǎn)力產(chǎn)生的殘差力。若殘差力足夠小，達(dá)到收斂標(biāo)準(zhǔn)，則該弧長步的迭代收斂，下一個弧長步的迭代開始。若達(dá)不到收斂標(biāo)準(zhǔn)，則繼續(xù)迭代求解。

(16)

(17)

可求得：

(18)

這樣就確定了Ai+1的位置，其對應(yīng)的狀態(tài)為

(19)

更新變量，返回到式(10)，準(zhǔn)備下一個迭代：

(20)

Δλi+1=Δλi+μ

(21)

i=i+1

(22)

(23)

5 算例

長度L=600 mm，半徑R=300 mm，厚度t=2 mm，彈性模量E=70 000 MPa，泊松比μ=0.3的圓筒殼，底部簡支，自由端施加軸向壓力，分析軸壓穩(wěn)定性。

5.1 根據(jù)工程計算公式求解

航天飛行器結(jié)構(gòu)設(shè)計應(yīng)用中，對大量圓柱殼體的軸壓穩(wěn)定性試驗(yàn)的數(shù)據(jù)進(jìn)行統(tǒng)計，歸納出如圖4的設(shè)計曲線，即軸壓臨界應(yīng)力系數(shù)k0～R/t曲線。根據(jù)圖4所示的曲線查取相應(yīng)的k0，再根據(jù)公式σcr=k0Et/R進(jìn)行計算。

由圖4可以看出，根據(jù)設(shè)計曲線求得的臨界軸壓非常保守。

圖4 軸壓臨界應(yīng)力系數(shù)k0～R/t曲線

5.2 有限元法

圖5 特征值屈曲模態(tài)

由計算分析看出，非線性屈曲的臨界應(yīng)力T?明顯低于特征值屈曲求得的臨界應(yīng)力T″，略高于經(jīng)驗(yàn)公式求得的臨界應(yīng)力T′。非線性屈曲的臨界應(yīng)力T?和經(jīng)驗(yàn)公式的臨界應(yīng)力T′之比為T?/T′=1.1/1，二者比較接近。

改變圓筒殼的長度和厚度，同樣可以進(jìn)行非線性屈曲的求解。幾種工況下有限元非線性屈曲分析求得的臨界軸壓和工程公式計算得到的臨界軸壓之比T?/T′如表1所示。由表1可見，有限元非線性屈曲的臨界軸壓計算結(jié)果和采用工程公式的計算結(jié)果吻合度較好。其中，t=3 mm時k0可取0.26。

圖6 非線性屈曲失穩(wěn)變形

圖7 載荷比例因子和弧長的關(guān)系

L=600 mmL=900 mmL=1 200 mmt=2 mm1.1/11.1/11.05/1t=3 mm1.01/11.04/11.02/1

對應(yīng)上述6種工況的線性特征值屈曲模態(tài)和非線性屈曲失穩(wěn)變形圖如圖8～圖13所示。每個圖的左圖為線性屈曲模態(tài)示意圖，云圖為位移，右圖為非線性屈曲失穩(wěn)示意圖，云圖為MISES應(yīng)力。

圖8 L=600 mm，t=2 mm時屈曲變形

圖9 L=600 mm，t=3 mm時屈曲變形

圖10 L=900 mm，t=2 mm時屈曲變形

圖11 L=900 mm，t=3 mm時屈曲變形

圖12 L=1 200 mm，t=2 mm時屈曲變形

圖13 L=1 200 mm，t=3 mm時屈曲變形

6 結(jié)論

1) 采用特征值屈曲模態(tài)缺陷和修正的Riks算法，對薄壁圓筒進(jìn)行了屈曲分析，并將分析結(jié)果和工程計算公式的分析結(jié)果進(jìn)行了對比，二者具有一致性。

2) 有限元法對任意結(jié)構(gòu)均能求解。對于無法用公式求得臨界軸壓的結(jié)構(gòu)，采用基于特征值屈曲模態(tài)缺陷的弧長法進(jìn)行非線性屈曲求解，可以得到合理的結(jié)果。

3) 采用特征值屈曲模態(tài)缺陷進(jìn)行分析，不必事先知道缺陷的具體位置和缺陷的性質(zhì)，在特征值屈曲的基礎(chǔ)上引入厚度的10%作為缺陷參數(shù)，采用弧長法進(jìn)行非線性屈曲分析，可以滿足工程設(shè)計需要，尤其是在工程研制的方案階段，在沒有真實(shí)產(chǎn)品用來進(jìn)行試驗(yàn)的情況下，采用這種方法進(jìn)行前瞻性的設(shè)計和分析，具有重要的意義。