王驍峰，段 毅，袁銳之

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

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

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

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

1 非線性屈曲的描述

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

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

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

2 特征值屈曲

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

(1)

(2)

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

(3)

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

3 廣義弧長(zhǎng)法求解非線性屈曲

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

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

(4)

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

4 修正的Riks算法

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

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

圖2 弧長(zhǎng)法求解非線性屈曲

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

(5)

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

(6)

(7)

從而

(8)

此處

(9)

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

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

Δλi=Δλ0

(10)

(11)

4.2 迭代過(guò)程

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

(12)

(13)

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

(14)

檢查迭代平衡

(15)

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

(16)

(17)

可求得：

(18)

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

(19)

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

(20)

Δλi+1=Δλi+μ

(21)

i=i+1

(22)

(23)

5 算例

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

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

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

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

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

5.2 有限元法

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

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

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

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

圖7 載荷比例因子和弧長(zhǎng)的關(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

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

圖8 L=600 mm，t=2 mm時(shí)屈曲變形

圖9 L=600 mm，t=3 mm時(shí)屈曲變形

圖10 L=900 mm，t=2 mm時(shí)屈曲變形

圖11 L=900 mm，t=3 mm時(shí)屈曲變形

圖12 L=1 200 mm，t=2 mm時(shí)屈曲變形

圖13 L=1 200 mm，t=3 mm時(shí)屈曲變形

6 結(jié)論

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

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

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