基于弧長法的南瓜型超壓氣球非線性后屈曲分析

卜亞樓， 楊燕初*， 蔡榕， 趙榮， 張向強

(1.中國科學院大學航空宇航學院， 北京 100040； 2.中國科學院空天信息創新研究院， 北京 100094)

南瓜型超壓氣球是目前長航時浮空飛行平臺的主流設計形式。其通過在歐拉體球形[1]表面均勻鋪設加強筋、使球膜形成局部高曲率鼓包的方式來減小球膜應力，提高承壓能力。因形狀似南瓜，因此被稱作南瓜型超壓氣球。南瓜球具有駐空高度穩定、飛行時間長、效費比高的特點，是執行臨近空間科學試驗和高空飛行任務最穩定的平臺之一[2]。南瓜球雖然在承壓能力和應力分布方面具有明顯優勢，但在某些情況下卻呈現出扭曲、不對稱的構型，使得南瓜球的穩定性研究和分析成為其結構設計的重要問題。

針對南瓜球的失穩現象，Calladine[3]最早從幾何穩定性的角度開展研究。通過將球膜的幾何剛度與歐拉梁的抗彎剛度進行類比，得到球膜幅數n與凸出角α(球膜局部凸出鼓包的弧長所對應的圓心角)之間的穩定性判據。但在隨后的研究中發現，Calladine[3]的結論更適用于按照恒定凸出角度(constant bulge angle, CA)設計的氣球。Smith等[4]認為按照恒定凸出半徑(constant bulge radius, CR)設計的氣球展開穩定性更好；而Schur[5]通過地面試驗推斷，過量的幅寬可能導致南瓜球的展開不穩定性，并指出按照恒定凸出半徑設計的氣球更有利于展開；Baginski等[6-9]通過能量最小原理計算氣球的平衡構型。對成型氣球進行穩定性分析發現，如果氣球處于不穩定的平衡狀態，那么在正常上升過程中可能不會完全展開。

有限元技術的發展使得利用數值方法來求解結構的失穩問題成為可能。文獻[10-11]使用基于動態松弛求解方法的內部有限元程序套件inTENS來模擬南瓜球的失穩行為。該方法基于顯式求解技術，其時間步動態分析結合了自動動能阻尼控制技術，通過給模型設定內部壓力和相應的薄膜應力進行數值預測。文獻[10-11]使用該程序套件對直徑為10 m的南瓜球進行求解，得到了內部壓力為690 Pa時的扭曲形態，數值預測結果與試驗相比有良好的相關性；Deng[12]基于ABAQUS/Explicit計算了部分充氣氣球的平衡構型。其中考慮了球膜材料的非線性行為，對正交各向異性黏彈性材料進行了分析和建模，并基于變泊松比的起皺準則模擬了球膜的褶皺行為，獲得的仿真結果與試驗相比有較好的一致性。Deng[12]還在仿真的基礎上結合量綱分析構建了南瓜球穩定性判據的經驗公式，是南瓜球穩定性研究上的一大進展；Pagitz等[13]認為南瓜球扭曲的穩定平衡構型是一個涉及平衡路徑分叉屈曲的問題，他們通過求解整體切線剛度矩陣的特征值和特征向量來研究失穩壓強和失穩形態；Xu等[14-15]利用ABAQUS軟件開展特征值屈曲分析，預測了氣球的臨界失穩壓強和相應的屈曲模態，對可能影響南瓜球穩定性的因素進行靈敏度分析，獲得了重要的研究結論；之后在特征值屈曲分析的基礎上基于*static, stabilize開展非線性后屈曲分析以模擬氣球超過臨界壓力后的響應。

Koiter[16]提出了考慮初始缺陷的非完善結構穩定性的一般準則和缺陷敏感度的概念，指出實際結構不可避免地存在初始缺陷，這會降低結構的屈曲載荷。具體地，Baginski等[7]指出南瓜球制造的不精確性、球膜和加強筋材料的不均勻性，以及外界環境的變化(如溫度)對穩定性有很大影響。這表明，初始缺陷對南瓜球穩定性的影響不可忽略。Xu等[15]中雖構造了初始缺陷開展非線性后屈曲分析，但其側重點在于研究氣球超過臨界壓力后的響應，無法獲得缺陷模式、缺陷幅值等對南瓜球極限承壓能力的影響。因此，目前對于初始缺陷如何影響南瓜球穩定性的問題僅停留在定性認識上，尚未開展深入研究。

較多的薄膜結構分析采用顯式準靜態算法[17]，而南瓜球由于整體結構尺寸較大、球膜具有較強的幾何非線性且球膜和加強筋彈性模量相差懸殊，傳統計算方法面臨著計算量大、計算不易收斂的問題。更重要的是，顯式準靜態算法無法獲得缺陷模式、缺陷幅值等對南瓜球極限承壓能力的影響。而基于弧長法的非線性后屈曲分析是研究結構缺陷敏感性的重要方法，其在追蹤失穩路徑、捕捉極限承壓力方面具有明顯優勢。文獻[18]采用弧長法對膜結構雷達罩進行了穩定性研究，取得了較好的結果，是弧長法在膜結構領域的最新應用。此外，該方法在球殼受外壓屈曲[19-20]、網殼結構屈曲[21-22]、復合材料加筋壁板受軸壓屈曲[23-25]等領域也有廣泛應用。

因此，在前人研究的基礎上，現將弧長法應用到南瓜球穩定性研究上。基于該方法既能獲得南瓜球的缺陷敏感特性，也能獲得南瓜球前后屈曲的全部信息,是南瓜球穩定性研究的一種新的方法。采用弧長法對受內壓南瓜球的非線性屈曲過程進行數值模擬，基于特征值屈曲模態位移構建了兩種缺陷模式，對模態缺陷下受壓南瓜球的穩定性進行研究。通過對比不同缺陷幅值和缺陷模式下南瓜球的臨界載荷，研究南瓜球對幾何缺陷的敏感性。針對弧長法下降段難以收斂的問題，使用重啟動方法結合顯式動力學計算南瓜球屈曲后的響應，獲得南瓜球的失穩形態和應力分布，為南瓜球在工程領域的設計制造提供參考。

1 穩定性分析的有限元法

1.1 線性特征值屈曲分析

特征值屈曲分析是一種線性攝動分析，它基于小位移線彈性理論，可以計算包含桁架單元、梁單元、板殼單元、實體單元的結構臨界載荷系數和相應的屈曲模態。結構的靜力平衡方程為

K+λKGu=Q

(1)

式(1)中：K為結構的彈性剛度矩陣；KG為結構的幾何剛度矩陣；λ為特征值；u為結構的模態特征向量；Q為作用在結構上的載荷。

平衡方程[式(1)]失穩的條件是方程存在奇異性，即等效剛度矩陣的行列式的值為零。

|K+λKG|=0

(2)

線性屈曲分析即為求解式(2)的特征值，所求得的特征值即為臨界載荷系數。將某一階屈曲模態特征值與施加的計算載荷相乘，就得到該階模態下結構的特征值屈曲臨界載荷。u為結構的模態特征向量，即屈曲模態位移，它預測了結構可能的失效形式。由于線性特征值屈曲以小位移小應變的線彈性理論為基礎，沒有考慮結構在受載過程中結構構型的變化，在外力施加的各個階段，平衡方程始終建立在初始構型上，忽略了實際加載過程中結構剛度矩陣的變化。因此，線性特征值屈曲分析計算結果偏保守，對于缺陷敏感結構需要進行非線性后屈曲分析。

1.2 弧長法

弧長法是結構非線性分析中計算穩定高效的一種數值方法，能夠有效地分析結構非線性前后屈曲并跟蹤后屈曲路徑，從而獲得結構失穩前后的全部信息，并有效捕捉極限載荷。

在非線性后屈曲分析中，針對Newton-Raphson迭代方法在屈曲臨界點附近結構剛度接近奇異時迭代不易收斂的情況，通過結合修正的弧長法可以穩定求解通過極值點的結構分析全過程。假設載荷成比例加載，迭代方程為

(3)

將位移增量ΔUi與載荷比例因子λt+Δt通過弧長半徑L聯系起來可得到求解λt+Δt的約束方程，即

ω(Δλi)2‖PW‖2+β‖ΔUi‖2=(ΔL)2

(4)

式(4)中：‖PW‖、‖ΔUi‖為外載荷向量和位移增量的二范數；Δλi為第i次迭代的載荷比例因子增量；ΔL為弧長增量的半徑；ω、β為尺度因子，在不同的弧長控制方法中取值不同。

2 有限元建模與特征值屈曲分析

2.1 數值模型建模

取某一南瓜球作為算例，其幾何參數如表1所示。

表1 南瓜球幾何參數

圖1為赤道平面處相鄰加強筋之間球膜的截面形狀。氣球的加強筋位于相鄰兩幅球膜的焊縫處，氣球的赤道半徑用R表示。對于n幅球膜的情況，赤道對角為θ=2π/n。

s為鼓包的弧長；c為弦長；r為凸出半徑；α為凸出角圖1 赤道平面處的凸出截面形狀Fig.1 Lobe geometry at equator

根據幾何關系，弦長c可用赤道半徑和赤道對角表示：

(5)

此外，弦長也可用凸出半徑和凸出角表示：

(6)

根據式(5)和式(6)，可得出赤道半徑與凸出半徑的關系為

(7)

鼓包弧長為

s=αr

(8)

由于鼓包的圓弧與凸出角或凸出半徑有關，目前南瓜球的設計主要采用兩種方式，一是凸出半徑沿子午線保持恒定的凸出半徑設計，二是凸出角沿子午線保持恒定的凸出角度設計。祝榕辰[26]證明了恒定凸出半徑的設計方法能使球膜表面的應力分布更加均勻；文獻[5]從展開穩定性的角度論證了恒定凸出半徑設計的氣球具有更好的穩定性。因此，恒定凸出半徑設計方法得到了較廣泛應用。本文的計算模型采用恒定凸出半徑設計。根據以上幾何參數及幾何關系完成建模，三維模型圖如圖2所示。

圖2 南瓜球三維模型示意圖Fig.2 Schematic layout of a 3D pumpkin balloon

氣球球膜和加強筋對其所用的材料有多種特性要求，隨著材料科學的發展，球膜和加強筋的特性在不斷改善。線性低密度聚乙烯以其優越的性質取代了以前的低密度聚乙烯作為球膜材料；PBO(p-phenylene benzobisoxazole)以其優越的力學性能、熱穩定性和優異的強度重量比被選擇作為加強筋的材料。目前，世界上幾個開展高空氣球項目的國家均在采用這種材料，如ULDB(ultra long duration balloon)項目。具體材料屬性如表2所示。

南瓜球的線彈性各向同性球膜采用M3D4R膜單元建模。這是由于膜單元僅有平面內剛度，沒有彎曲剛度，且M3D4R是一個4節點四邊形單元，為減縮積分單元，具有計算高效的優點；桁架單元通常用于模擬僅在軸向承受拉伸或壓縮應力的細長結構，氣球中的加強筋雖然無法承受壓應力，但在超壓狀態下所承受的拉應力與桁架單元沿軸向傳遞的拉伸特性是相同的，因此采用桁架單元T3D2對加強筋進行建模，其中加強筋的剛度定義為彈性模量和橫截面的乘積。

南瓜球超壓時，載荷會從球膜轉移到位于焊縫處的加強筋上，并連接到位于氣球頂部和底部的剛性法蘭盤上。位于焊縫處的加強筋通過封套約束在球膜表面，加強筋和球膜之間沿軸向為滑動摩擦約束。但實際數值仿真中，由于筋膜之間的滑動摩擦約束關系不易建模，且筋膜之間容易相互穿透導致計算不收斂，因此通常將筋膜建模為共節點約束以代替滑動摩擦約束來傳遞軸向力。趙榮等[27]已通過數值仿真證明了加強筋和膜單元之間采用共節點約束與采用滑動摩擦約束對球膜承力并無很大差別；文獻[28]同樣采用共節點約束開展靜力學分析，取得了較為合理的結果。因此，為了建模簡單起見，本研究同樣采用共節點約束。頂部法蘭盤和底部法蘭盤均為全約束。此外，為了獲得較好的網格質量，本分析中沒有構建法蘭盤部件，而是以桁架單元代替法蘭，這樣做可以提高計算的收斂性而不影響計算結果。單元數量對于模擬結果的準確度影響很大，需要進行無關性驗證，一般認為單元數量的選取應該使前后兩次分析結果相對浮動量小于等于1%。經計算比較，最后選擇的劃分單元數量為T3D2：13 600個，M3D4R：64 600個。取單位載荷1 Pa為計算載荷，以均布力形式施加在氣球內表面。

表2 材料屬性

2.2 線性特征值屈曲分析

首先對南瓜球進行線性特征值屈曲分析，得到屈曲壓強和模態位移。在非線性后屈曲分析中，需要用模態位移構造初始幾何缺陷。

計算得到前8階屈曲壓強和屈曲模態，結果如表3所示。從表3中看出，在前8階屈曲模態中存在著模態形狀相似、特征壓強相近的4組相鄰模態，如1階和2階、3階和4階，其屈曲模態構型除了在位移幅值上有微小的差別外，相近特征值下的兩個模態為繞豎直軸旋轉約15°(圖3)，說明即使在相等的外載荷作用下，南瓜球仍存在多種可能的屈曲位移趨勢，證明南瓜球具有分支屈曲的特性。

圖4顯示了南瓜球前4組相近特征值對應的屈曲模態。需要注意的是屈曲模態u為歸一化向量，最大位移分量為1.0，這并不代表在臨界載荷下變形的實際大小，但能夠表明結構的可能失穩模式。屈曲模態構型表現出整體對稱的n上n下(n為自然數)形態。高階模態對應更大的n值；隨著n的增大，需要更大的臨界壓強通過分叉加載點，但臨界壓強與n值并非簡單的線性關系。n不超過4時，相鄰特征壓強相差較小；當n大于4時，相鄰特征壓強表現出顯著的差距。

表3 南瓜球特征值屈曲分析結果

圖3 具有相近特征值的兩個屈曲模態對比圖Fig.3 Comparison diagrams of two buckling modes with similar eigenvalues

圖4 前四組特征模態的總位移等值線圖Fig.4 Contourplots of total displacements of the first four eigenmodes

3 非線性后屈曲分析

3.1 初始缺陷的建立

氣球受壓時，在屈曲前主要受均勻的面內力，保持其球形均勻擴張的平衡狀態；到達屈曲臨界點時，有多個可能的失穩分支路徑；當受到擾動時，將沿某一路徑發生變形失穩，進入到新的平衡狀態，屬于分支點屈曲問題。基于線性屈曲分析的模態位移構造初始缺陷分布用以誘導結構的后屈曲行為，利用弧長法進行非線性后屈曲分析，獲得南瓜球的缺陷敏感特性。

其中應用式(9)來表述缺陷幾何的擾動網格：

(9)

式(9)中：ψi為第i階模態形狀；θi為相應的縮放因子，即缺陷比例系數，表示對模態位移進行縮放；Δxi為縮放后的模態位移。

由式(9)可知，構建初始缺陷的方式是多樣的，既可以對某一階模態位移進行縮放，也可以對多階模態位移縮放后進行線性組合。因此在進行非線性后屈曲分析時，采用了兩種方法來構造初始缺陷。第一種是臨界缺陷模態法，即基于4上4下臨界模態位移的缺陷氣球。該理論認為，臨界屈曲模態對應最小勢能狀態，當結構缺陷分布形式與臨界模態吻合時，結構最容易發生屈曲，即最低階模態缺陷是最不利的幾何缺陷構型。運用臨界缺陷模態法對模型進行分析時，先后引入5種缺陷幅值進行對比分析，最小缺陷比例系數為0.01。

圖5顯示了帶有1%初始缺陷的氣球的幾何形狀，可以看出，氣球輕微變形，但足以觸發后屈曲行為。由于缺陷均是基于臨界模態位移，所以5種缺陷幅值所對應的構型一致，只是初始撓度不同。

圖5 基于臨界屈曲模態的初始缺陷氣球Fig.5 Balloon with initial imperfection based on first buckling eigenmode

另一種是組合缺陷模態法。Triantafyllidis等[29]發現，在一些不穩定結構中，由于其特殊的幾何形狀，在相同的載荷水平下會出現整體屈曲和局部屈曲。對于這些結構，不穩定點對初始缺陷的形狀非常敏感。該理論認為，不止一階屈曲模態用來構造初始缺陷，多階屈曲模態的組合也可用來構造幾何缺陷。因此，本文的第二種初始缺陷的構造方法是使用前8階模態位移進行線性組合來觸發后屈曲行為，其中每種模態的缺陷比例系數均為0.01。

3.2 缺陷敏感度分析

采用弧長法對缺陷氣球進行非線性有限元分析。圖6為某一缺陷幅值下載荷和應變能隨迭代弧長變化的曲線。當結構失穩時，載荷-位移曲線出現負剛度，對應的載荷-弧長曲線出現下降段，結構同時釋放應變能維持平衡。由圖6可知，載荷和應變能在同一時刻發生轉折，說明結構發生失穩。

圖6 載荷和應變能隨迭代弧長變化曲線Fig.6 Curve of load and strain energy with iteration arc length

圖7和圖8分別為臨界缺陷模態法和組合缺陷模態法在5種缺陷幅值下氣球的載荷-弧長曲線。屈曲載荷結果對比如表4所示。

圖7 基于臨界缺陷模態法不同缺陷幅值下的載荷- 弧長曲線Fig.7 Load-arc length curves under different defect amplitudes based on critical defect mode method

圖8 基于組合缺陷模態法不同缺陷幅值下的載荷- 弧長曲線Fig.8 Load-arc length curves under different defect amplitudes based on combined defect mode method

由表4可知，缺陷狀態下非線性屈曲臨界載荷小于線性特征值屈曲臨界載荷，特征值屈曲分析結果偏保守。因此，工程實際中有必要對南瓜球進行非線性后屈曲分析。

表4 屈曲載荷結果對比

對比相同缺陷模式下不同缺陷幅值所對應的屈曲載荷發現，臨界載荷隨缺陷幅值的增加而減小。假設實際制造的南瓜球其初始缺陷符合臨界屈曲模態構型1%的情況，其屈曲臨界載荷為224.5 Pa。當缺陷比例系數為0.02時，其屈曲臨界載荷為215 Pa，承載力相對下降約4.2%。當缺陷比例系數為0.05時，承載力相對下降13%，說明幾何缺陷幅值對屈曲臨界載荷有顯著影響。

圖9 充氣過程中位移等高線圖Fig.9 Contours of displacement during inflation

此外，缺陷幅值相同時，不同缺陷模式對應的非線性臨界屈曲載荷各不相同。基于組合特征缺陷模態法構造的幾何缺陷所對應的屈曲載荷要低于同幅值下的臨界缺陷模態法。同樣地，假設實際制造的南瓜球其初始缺陷符合組合屈曲模態構型1%的情況，其屈曲臨界載荷為212 Pa。當缺陷比例系數為0.02時，其屈曲臨界載荷為193.5 Pa，承載力相對下降約8.7%，當缺陷比例系數為0.05時，承載力相對下降53.5%。下降幅值高于相同缺陷比例系數下的臨界缺陷模態情況，且隨著比例系數的增大，結構承壓力大幅度下降，這表明：①缺陷幅值和缺陷模式對失穩載荷有顯著影響，南瓜球是一種缺陷敏感性結構；②第一階模態缺陷不一定是最不利的幾何缺陷構型，求解的臨界載荷有可能高估結構的失穩載荷；組合模態缺陷也有出現的可能，其對應的屈曲臨界載荷有可能更低，對南瓜球屈曲特性的研究應考慮多種模態缺陷情況，以便找出最低承載力對應的缺陷模式。

3.3 失穩形態和應力分布分析

針對弧長法下降段難以收斂的情況，利用重啟動方法結合顯示動力學計算了兩種缺陷模式下屈曲之后的失穩形態和應力分布隨載荷的變化情況。需要注意的是，由于球膜的大變形，球膜之間會發生接觸，需對接觸進行設置。

對于單個臨界缺陷模態的初始構型，發現氣球最終穩定為4上4下的整體屈曲構型。圖9顯示了氣球從初始缺陷構型充氣到500 Pa的過程中位移變化情況。在充氣剛開始時(100 Pa之前)，氣球僅發生輕微變形，位移變化較小，位移云圖開始出現不均勻性；繼續充氣，壓力從100 Pa變化到230 Pa的過程中氣球仍未出現明顯變形；從230 Pa開始，變形速率顯著增大；在300 Pa之后再次減小。氣球變形速率變化的轉折點位于230 Pa附近，大約為非線性后屈曲的臨界失穩壓強。在整個充氣過程中，氣球的缺陷模式被逐漸放大。

球膜的Mises應力云圖如圖10所示。高應力區位于氣球凸起的波峰，這可能是試驗氣球在高壓下爆裂的區域。低應力區位于波峰之間的接觸區域。觀察應力圖例條可知，應力值出現兩極分化，應力分布極不均勻，且最大應力值(近似取22.1 MPa)遠高于未失穩狀態時的最大應力值(約 4 MPa，可參考文獻[27])。南瓜球在承壓能力、應力分布方面的優勢，在面對失穩問題時消失。

圖10 Mises應力云圖Fig.10 Mises stress contour

對于組合缺陷模態的構型，氣球在500 Pa壓強下呈現局部變形，如圖11所示。與整體變形構型(4上4下)不同，圖11顯示了一種非均勻變形模式，其中幾幅球膜形成局部塌陷，而球體表面的其他部分比較均勻。與整體變形模式中應力分布不同，高應力區位于局部塌陷的區域，并且應力水平高于整體變形構型中的應力水平。

中外南瓜球的地面試驗曾多次出現以上兩種失穩構型。一種是整體變形模式，氣球具有循環對稱的變形形狀，如4上4下；另一種是局部變形模式，美國國家航空航天局(National Aeronautics and Space Administration，NASA)在氣球試驗中多次發生S形裂縫的失穩形態，中國科學院在地面試驗中曾出現過局部塌陷(圖12)的失穩形態。兩種不同的失穩形態可能是由于初始缺陷模式不同造成的。

圖11 Mises應力云圖Fig.11 Mises stress contour

圖12 地面試驗中南瓜球出現局部失穩Fig.12 Local instability of pumpkin balloon in ground test

在整體屈曲和局部屈曲構型中，加強筋出現嚴重變形。圖13顯示了裂縫區域周圍嚴重變形的加強筋形狀。其應力云圖顯示，無論是整體屈曲還是局部屈曲，在球膜接觸區域內，加強筋應力要相對偏低，并且其他區域內的加強筋具有幾乎均勻的應力值。這可以解釋為，在球膜的接觸區域，加強筋均位于較低位置，處于相對“松弛”狀態，這將導致應力降低。

圖13 加強筋應力分布Fig.13 Stress distribution of tendons

4 結論

通過有限元中的弧長法對基于屈曲模態位移構造的缺陷南瓜球進行非線性后屈曲分析，得到以下結論。

(1)南瓜球是缺陷敏感性結構，其屈曲特性的數值分析需要考慮初始缺陷。南瓜球的屈曲臨界載荷隨缺陷幅值的增大而降低，對于臨界缺陷模態法，在缺陷系數為0.02時，南瓜球的承壓力相對理想氣球下降約4.2%，表明提高南瓜型超壓氣球穩定性需要提高氣球的加工精度。

(2)采用臨界缺陷模態法和組合缺陷模態法模擬氣球初始缺陷構型，相同缺陷幅值下組合缺陷模態法比臨界缺陷模態法求得的屈曲臨界載荷低了至少4.5%，表明第一階模態缺陷不一定是最不利的幾何缺陷構型，對南瓜球屈曲特性的研究應考慮多種缺陷模式，避免單一缺陷模式求解的承壓力高估結構的穩定性。

(3)不同的缺陷模式在相同的壓強下形成的失穩構型不同。由臨界缺陷模態法模擬的初始缺陷構型，在500 Pa時穩定為整體失穩；由組合缺陷模態法模擬的初始缺陷構型，在500 Pa時穩定為局部失穩。

(4)南瓜球在內壓下失穩后應力分布極不均勻。整體屈曲構型中高應力區位于凸起的波峰上，局部屈曲構型中高應力區位于局部失穩區域，且失穩狀態下高應力區的應力水平遠高于相同載荷下未失穩狀態的應力水平(前者約為后者的5倍)。失穩現象大大降低了氣球的承壓能力，是南瓜球設計過程中急需解決的問題。