薛鵬聰，劉 鋮,2,3，水小平

（1. 北京理工大學 宇航學院，北京 100081；2. 中物院高性能數(shù)值模擬軟件中心，北京 100088；3. 北京應用物理與計算數(shù)學研究所，北京 100088）

引 言

膜結構能夠在保證結構可靠性的同時耗費更少的材料包圍更大的空間，并具有自重更輕、造型優(yōu)美、便于拆卸等優(yōu)點，因此在航天工程領域得到了廣泛的關注，被期望作為執(zhí)行深空探測任務的新一代創(chuàng)新空間結構，如薄膜太陽帆推進系統(tǒng)、通訊衛(wèi)星中大型薄膜天線反射器及太陽能電站中電池陣列結構等[1]。這類空間結構在發(fā)射前通常處于折疊狀態(tài)，入軌后在驅動裝置的作用下展開，并進入工作狀態(tài)。然而對于在太空中展開太陽帆這類具有高收納比的大型薄膜結構，無疑是一項復雜的操作，且展開后的有效面積及膜面精度均會對其所承擔的功能產生一定影響。因此，進入太空中的薄膜航天器能否順利展開至理想工作狀態(tài)是這類航天器折展系統(tǒng)設計中的關鍵技術?？紤]到地面實驗的成本過大，且難以針對大型薄膜結構提供在軌微重力真空環(huán)境[2]。于是，采用高效的數(shù)值計算方法對空間帆膜結構的動力學展開過程進行仿真模擬是此類薄膜航天器結構設計的重要手段。Shirasawa等[3]采用多質點模型對“伊卡洛斯”（Icarus）太陽帆展開動力學進行數(shù)值模擬，采用由彈簧和阻尼器連接的質點代替?zhèn)鹘y(tǒng)薄膜單元，近似簡化了膜結構建模，降低了計算成本。Miyazaki[4]提出了一種剛度退化的褶皺模型，通過引入剛度折減系數(shù)描述薄膜壓縮時的本構關系，建立了薄膜太陽帆的柔性多體系統(tǒng)有限元模型，并將其應用于薄膜太陽帆的展開動力學分析。

針對類似薄膜結構的柔性多體系統(tǒng)的大轉動、大變形問題，柔性多體系統(tǒng)動力學領域提出了絕對節(jié)點坐標方法（Absolute Nodal Coordinate Formulation，ANCF）[5]。基于ANCF方法，國內學者針對薄膜太陽帆等多柔體航天器的展開動力學進行了較為深入的研究[6-9]。趙將等[7]采用ANCF薄板單元，對考慮粘彈性阻尼的“伊卡洛斯”薄膜太陽帆進行了展開動力學分析。Liu等[8]系統(tǒng)地對ANCF板/殼單元進行了研究，結合SRM模型提出了縮減的ANCF薄膜單元，并通過薄膜太陽帆的自旋展開動力學模擬驗證了所提出方法的有效性。然而，由于薄膜結構在宏觀響應中拉壓模量呈現(xiàn)的巨大差異，系統(tǒng)非線性特性顯著，導致薄膜結構在靜力學/動力學模擬中極易發(fā)生計算不穩(wěn)定現(xiàn)象。

針對拉壓不同模量彈性結構，前蘇聯(lián)學者阿姆巴book=218,ebook=122爾楚米揚（Ambartsumyan）[10]系統(tǒng)地提出了不同模量彈性理論，并給出了其本構關系的基本假定，推導了二維平面問題和三維空間問題的廣義彈性定律。之后許多研究者基于不同角度對Ambartsumyan的不同模量彈性理論進行了發(fā)展和修正[11]。近年來，針對雙模量結構的數(shù)值算法得到了發(fā)展，張亮等[12]基于參數(shù)變分原理和有限元方法，在主應力方向上建立了統(tǒng)一的含參數(shù)變量的本構方程，將不同模量問題演化為一個易于求解的標準互補二次規(guī)劃問題，避免了在有限元迭代過程中不斷根據應力狀態(tài)更新剛度矩陣，并將其應用于二維薄膜的起褶分析，呈現(xiàn)出不錯的收斂效率。杜宗亮等[13]基于不同模量彈性理論推導了一種具有漸進二次收斂效率的切線剛度矩陣，通過將其應用于二維薄膜的應力水平和起褶區(qū)域的預測，驗證了該方法在處理不同模量薄膜結構時的高效性、穩(wěn)定性。

本文旨在基于ANCF建模方法，整合張力場與不同模量彈性理論[10]，提出一種考慮不同模量的ANCF薄膜單元，以實現(xiàn)空間薄膜結構的高效動力學模擬。為提高單元在處理拉壓不同模量薄膜問題時的算法穩(wěn)定性，推導了單元精確的切線剛度矩陣，并通過若干靜力學算例驗證了單元的高效性。此外，針對一正六邊形薄膜太陽帆，建立其展開過程多柔體系統(tǒng)動力學模型，分析了薄膜太陽帆的結構參數(shù)對其展開過程的影響，數(shù)值仿真結果為這類空間帆膜結構的設計和優(yōu)化提供了理論指導。

1 基于ANCF描述的薄膜單元

如圖1所示，點P(ξ,η)為初始構形中薄膜單元上一點，該點的位置矢量r0(ξ,η)表示為

book=219,ebook=123

2 基于不同模量彈性理論的彈性力及切線剛度

2.1 薄膜張力場理論與不同模量彈性理論

張力場理論[16]認為薄膜的彎曲剛度為零，當薄膜內出現(xiàn)壓應力時，一般通過面外變形即薄膜起褶的形式釋放壓應力，張力場理論將其處理為薄膜的面內收縮。起褶后的薄膜處于單軸拉伸應力狀態(tài)，褶皺的方向平行于大主應力的方向，且垂直于褶皺方向的小主應力為零。

book=220,ebook=124

2.2 考慮不同模量的ANCF薄膜單元的切線剛度

3 基于廣義-α算法的多體系統(tǒng)動力學方程

圖2給出了薄膜結構的動力學分析流程圖，在預處理模塊中，根據薄膜結構的初始構形以及材料參數(shù)和單元信息等建立了精確的ANCF薄膜單元，給出了系統(tǒng)的常數(shù)質量矩陣。然后，采用高斯消元法消除冗余book=221,ebook=125約束，得到系統(tǒng)的稀疏矩陣。通過假定初始時刻系統(tǒng)的應力場確定初始彈性矩陣，并組裝系統(tǒng)的彈性力和雅克比矩陣，然后代入到式（36）～（39）的廣義-α算法的求解過程中，根據求解結果更新系統(tǒng)位移場和彈性矩陣，如此循環(huán)迭代直至達到給定的收斂準則，求解完成。

圖2 膜結構動力學分析計算方法Fig. 2 Dynamic analysis and calculation method of membrane structure

4 數(shù)值仿真驗證

為便于區(qū)分，將式（29）對應的線性近似的切線剛度模型記為ANCF LM，將式（30）對應的切線剛度模型記為ANCF TM?；诒疚奶岢龅挠嬎惴椒?，本節(jié)主要針對薄膜結構的一些經典案例進行了考察。4.1節(jié)、4.2節(jié)和4.3節(jié)分別驗證文所提出的計算框架在求解不同模量小變形、大變形，以及面外變形問題時的高效性和精確性。

4.1 基于ANCF的薄膜單元小變形分析

圖3所示為一個驗證薄膜起褶模型各種數(shù)值算法的經典算例，矩形薄膜寬度H，厚度t，其上下邊界施加有大小σ0的均勻預應力，左右兩端施加有軸向集中力P=σ0tH，同時左右邊界均施加有集中力矩M。隨著力矩的增加，膜的下邊緣開始出現(xiàn)一定寬度的褶皺區(qū)域?？紤]采用ANCF DM薄膜單元計算框架，對上述薄膜結構進行了數(shù)值模擬，仿真過程中，不考慮薄膜的壓縮模量。

圖3 純彎曲矩形薄膜受力示意圖Fig. 3 Stress diagram of Pure Bending Rectangular membrane

文獻[19]針對上述模型，采用解析方法對其起褶區(qū)域以及面內變形進行了研究，并得出一系列關系式。其中，式（40）給出了矩形膜起褶區(qū)域的寬度b與外力矩M的解析關系式

式（41）給出了不同工況下，薄膜右邊界上一點的柯西應力與外力矩的解析關系式

式（42）給出了薄膜變形后曲率κ與外力矩M之間的解析關系式

圖4給出了起褶區(qū)域寬度的解析解和數(shù)值試驗結果，不難看出即使在起褶區(qū)域超過90%時，數(shù)值模擬結果仍然具有很高的精度。

圖4 褶皺帶寬與外力矩的關系Fig. 4 Relationship between fold bandwidth and external torque

book=222,ebook=126圖5則給出了不同外力矩作用下，右邊界上一點柯西應力的解析解與數(shù)值試驗結果，兩者的吻合性比較好。

圖5 右邊界水平應力與外力矩的關系Fig. 5 Relationship between horizontal stress at right boundary and external moment

圖6給出了薄膜結構曲率的解析解和數(shù)值解隨外力矩的變化趨勢，分析得知，兩者的趨勢是非常一致的。

圖6 純彎曲矩形薄膜的彎矩–曲率關系Fig. 6 Moment-curvature relationship of rectangular membrane

通過上述各項數(shù)值試驗結果與解析結果的對比可以得出，本文所提出的方法在研究薄膜起褶的問題時具有很高的精度。

4.2 基于ANCF的不同模量結構大變形分析

如圖7所示為一L形梁結構，其左端為固定約束，自由端受水平方向集中力F= 40 000 N。固定其受拉方向的彈性模量為Et= 3×107Pa，泊松比為vt= 0.3，依次分別取壓縮模量為Ec=3×107Pa、Ec=1.5×107Pa、Ec= 6 ×106Pa、Ec=3×106Pa 以及Ec= 1.5 × 106Pa 時的工況，對應的拉壓比分別為Et/Ec=1、2、5、10、20。

為驗證本文所提出計算方法的有效性，本節(jié)基于ABAQUS平臺的擴展程序接口進行了二次開發(fā)，并將UMAT的結果作為ANCF描述的不同模量問題的參考解。需要指出的是，UMAT這里采用CPS4I平面應力四邊形單元和自適應時間增量步進行計算，在拉壓比較大的情況下需要較多的加載步。而本文所提出的ANCF TM模型僅需一個增量步即可得到收斂解。

圖7 L形梁結構Fig. 7 L-shaped beam structure

圖8給出了不同拉壓模量比下，L形梁的變形后構形圖，表1給出了采用不同技術途徑在不同拉壓模量比條件下得到的節(jié)點P的水平位移以及計算收斂耗時。以上仿真結果表明，即使在求解拉壓模量相差20倍的大變形問題時，本文所提出的計算框架仍能高效收斂，且能夠保證較高的精度。

圖8 變形后構形圖Fig. 8 Configuration diagram after deformation

表1 節(jié)點P的水平位移與計算耗時Table 1 Horizontal displacement of nodePand calculation convergence time

book=223,ebook=127此外，如表2所示，在采用ANCF-LM模型求解不同模量大變形問題時，無法得到收斂的結果。值得注意的是，即使拉壓模量相同的條件下，ANCF-LM模型仍然無法得到收斂的結果，這是由于在相同模量時，該模型無法退化為經典的線性彈性本構模型。而本文提出的ANCF-TM模型求解過程中牛頓迭代收斂性良好，迭代步數(shù)和迭代耗時T隨拉壓模量比的增大略有增加。

表2 不同模型下節(jié)點P的求解效率Table 2 Solution efficiency of nodeP under different models

4.3 安全氣囊的充氣仿真過程

圖9所示為一方形安全氣囊的初始構形，薄膜材料的楊氏模量為588.0 Mpa，泊松比為0.4。方膜對角線長度為L= 1 200 mm，厚度t= 1mm?？紤]將氣囊的充氣過程近似為一個準靜態(tài)過程，氣囊在內部遞增的均勻壓力作用下發(fā)生變形，直至內部壓強P最終增至5 000 Pa，充氣過程結束。仿真過程中，不考慮薄膜的壓縮模量。

圖9 安全氣囊的初始構形示意圖Fig. 9 Schematic diagram of initial configuration of airbag

表3給出了基于本文所提出的計算方法在不同網格劃分密度下得到的安全氣囊上節(jié)點O、C、B的位移以及計算耗時T。結果表明，基于本文計算框架得到的位移收斂解與文獻[20]的研究結果基本吻合。此外，采用10×10的網格即可得到位移收斂解，相應的計算耗時隨著網格的加密有大幅增加。

表3 不同網格密度下的薄膜節(jié)點的位移Table 3 Displacement of membrane nodes under different mesh densitied

圖10出了安全氣囊充氣完成后的應力分布示意圖，根據應力準則[13]，對單元內每個積分點的應力狀態(tài)進行標記，用黑點標記的薄膜區(qū)域處于褶皺狀態(tài)，用圓標記的薄膜區(qū)域處于松弛狀態(tài)，其余未標記的區(qū)域則處于張緊狀態(tài)，這一結果與文獻[21]的結論基本一致。綜上結論再次驗證了本文所提出的計算框架在求解薄膜面外變形問題時的可行性。

圖10 安全氣囊變形后應力分布示意圖Fig. 10 Schematic diagram of stress distribution of airbag after inflation

5 薄膜太陽帆的自旋展開動力學模擬

作為一種新型空間可展開結構，能夠通過自旋展開并最終穩(wěn)定的大型帆膜結構在未來的空間探索技術應用中具有巨大潛力。能夠準確地對薄膜展開過程的動態(tài)響應進行預測和模擬是薄膜太陽帆展開系統(tǒng)設計中的關鍵技術之一。這里將采用本文所提出的計算框架對薄膜太陽帆的自旋展開過程進行分析，通過確定影響展開的關鍵因素來指導太陽帆的結構設計。

初始時刻，薄膜太陽帆按照圖11所示的折痕折疊后纏繞在與其同心的正六邊形彀輪之上，彀輪驅動后，太陽帆隨之旋轉，并在均布于其六個角上的集中質量塊的離心力作用下自旋展開，完全展開后的太陽帆是一個正六邊形。此時太陽帆的邊長為L= 0.635 m，厚度為t= 0.1 mm。薄膜材料的密度為ρ= 1 420 kg/m3，book=224,ebook=128楊氏模量為E= 2.5×109Pa，泊松比為v= 0.3。彀輪的邊長為d= 0.1 m，質量為m0=1.0 kg，相對于過質心的Z軸的慣性矩為Iz= 5.0×10–3kg m2。薄膜邊界角點處的集中質量塊的質量為m= 0.05 kg，初始時刻系統(tǒng)的角速度為w= 0.628 rad/s。

圖11 薄膜太陽帆的有限元模型Fig. 11 Finite element model of membrane solar sail

考慮到結構的對稱性，本文僅建立1/6模型，并施加旋轉邊界條件，該有限元模型中包括35個薄膜片，薄膜片之間通過旋轉鉸連接，具體建模方法可以參考文獻[8]。本文采用考慮拉壓不同模量的ANCF薄膜單元和第一類拉格朗日方程對太陽帆進行動力學精確建模，然后利用可控數(shù)值耗散的廣義-α算法對太陽帆的動力學展開進行求解分析。

5.1 壓縮剛度折減系數(shù)對展開過程的影響

在對太陽帆進行仿真時考察了3種不同的剛度折減的薄膜模型（Ec=s×Et），其中s為剛度折減系數(shù)。模型一為不考慮壓縮剛度的不可壓ANCF薄膜單元（s=0）；模型二為考慮微弱抗壓剛度的ANCF薄膜單元（s=10–5）；模型三同樣為引入微弱抗壓剛度的ANCF薄膜單元（s= 10–4）。

圖12給出了采用模型一進行仿真時薄膜太陽帆的展開過程圖，結果表明，此時的薄膜太陽帆無法完全展開并維持展開狀態(tài)，太陽帆在12 s左右展開至最大后又重新收縮至中心彀輪附近。

圖12 太陽帆的展開過程圖（s=0）Fig. 12 Deployment process diagram of solar sail （s=0）

圖13和圖14分別給出了采用模型一進行仿真時太陽帆展開過程中的能量與角動量演化過程。圖13表明，即使展開過程中薄膜的應變能出現(xiàn)高頻振蕩，但是系統(tǒng)的總能量仍然是守恒的。在12 s左右時太陽帆系統(tǒng)的應變能達到了峰值，此時太陽帆處于最大張緊狀態(tài)。圖14表明，相比于中心彀輪和集中質量塊，薄膜在整個仿真過程中的角動量僅發(fā)生小范圍的波動，在仿真時間達到12 s之后，太陽帆系統(tǒng)各部件的角動量均趨于穩(wěn)定。

圖13 太陽帆展開過程中系統(tǒng)的應變能及總能量（s=0）Fig. 13 Strain energy and total energy of the system during solar sail deployment （s= 0）

圖14 太陽帆展開過程中系統(tǒng)各部件的角動量（s=0）Fig. 14 Angular momentum of each component of the system during solar sail deployment （s=0）

圖15給出了采用模型二進行仿真時薄膜太陽帆的展開過程圖，此時的薄膜太陽帆同樣無法維持最大展開狀態(tài)，太陽帆在11 s左右展開到最大半徑，之后在邊緣集中質量區(qū)域出現(xiàn)較大幅度的面外變形。根據圖16book=225,ebook=129和圖17給出的太陽帆系統(tǒng)展開過程中的能量與和角動量的演化過程，可以得知，系統(tǒng)的應變能在12 s之后穩(wěn)定在一較小的范圍內。此外，在仿真時間達到11 s左右時，太陽帆系統(tǒng)各部件的角動量出現(xiàn)突變，之后逐漸穩(wěn)定，展開過程中系統(tǒng)的總能量和總角動量均是守恒的。

圖15 太陽帆的展開過程圖（s=10-5）Fig. 15 Deployment process diagram of solar sail（s=10-5）

圖16 太陽帆展開過程中系統(tǒng)的應變能及總能量（s=10-5）Fig. 16 Strain energy and total energy of the system during solar sail deployment（s=10-5）

圖17 太陽帆展開過程中系統(tǒng)各部件的角動量（s=10-5）Fig. 17 Angular momentum of each component of the system during solar sail deployment （s=10-5）

圖18給出了采用模型三進行仿真時薄膜太陽帆的展開過程圖，這時的薄膜太陽帆可以完全展開并維持展開狀態(tài)。太陽帆在13 s左右展開至最大半徑，之后跟隨彀輪一起穩(wěn)定旋轉，同時在邊緣集中質量區(qū)域發(fā)生微弱的面外變形。同時根據圖19和圖20分別給出其展開過程中的能量和角動量的演化過程可以得知，此時薄膜太陽帆系統(tǒng)的總能量和總角動量均是守恒的。

圖18 太陽帆的展開過程圖(s=10-4)Fig. 18 Deployment process diagram of solar sail(s=10-4)

圖19 太陽帆展開過程中系統(tǒng)的應變能及總能量（s=10-4）Fig. 19 Strain energy and total energy of the system during solar sail deployment（s=10-4）

圖20 太陽帆展開過程中系統(tǒng)各部件的角動量（s=10-4）Fig. 20 Angular momentum of each component of the system during solar sail deployment（s=10-4）

圖21和圖22給出了基于不同壓縮剛度折減系數(shù)的薄膜模型仿真得到的太陽帆系統(tǒng)各部件的動力學響應。圖21表明，當剛度折減系數(shù)為0和10–5時，太陽帆系統(tǒng)的彀輪和邊緣質量塊無法進入同步旋轉。而在剛度折減系數(shù)為10–4的情況下，在仿真時間達到13 s左右時彀輪和邊緣質量塊的旋轉角度開始相等，并進入同步旋轉。圖22表明，在考慮壓縮剛度的情況下，當剛book=226,ebook=130度折減系數(shù)較大時，薄膜太陽帆展開至一定半徑大小之后趨于穩(wěn)定。而當剛度折減系數(shù)較小或者不考慮壓縮剛度時，薄膜太陽帆無法維持最大展開狀態(tài)。綜上分析可知，在對薄膜太陽帆系統(tǒng)進行展開動力學分析時，考慮一定的壓縮剛度是必要的，這與文獻[4]得出的結果基本是一致的。

圖21 不同剛度折減系數(shù)時邊界質量塊和中心彀輪的轉角Fig. 21 Rotation angle of edge mass and central hub under different stiffness reduction factors

圖22 不同剛度折減系數(shù)時的太陽帆展開半徑Fig. 22 Deployment radius of solar sail under different stiffness reduction factors

book=227,ebook=131

5.2 集中質量和初始轉速對展開過程影響分析

為進一步確定薄膜太陽帆的一些結構參數(shù)對其展開動力學特性的影響，本文又分別仿真了當邊界集中質量和系統(tǒng)初始轉速變化時，太陽帆系統(tǒng)各部件的動力學響應。

圖23和圖24分別給出了太陽帆系統(tǒng)邊界集中質量為0.03、0.05、0.07 kg時其各部件的轉角演化過程和太陽帆展開半徑演化過程。分析得知，太陽帆分別在16、13、11 s左右展開至最大半徑，同時中心彀輪和邊界質量塊的旋轉角度開始趨于同步，展開過程的構形演化圖類似圖18。此外，綜合分析得知，邊界集中質量對太陽帆系統(tǒng)的最終展開半徑大小以及完全展開后發(fā)生的面外變形均會產生一定影響。

圖23 不同集中質量時邊界質量塊和中心彀輪的轉角Fig. 23 Rotation angle of edge mass and central hub under different concentrated masses

圖24 不同集中質量時的太陽帆展開半徑Fig. 24 Deployment radius of solar sail with different concentrated masses

圖25和圖26給出了系統(tǒng)初始轉速分別為0.471、0.628、0.785 rad/s時太陽帆展開過程中各部件的動力學響應。結果表明，太陽帆分別在16.5、13、10 s左右展開至最大半徑，同時中心彀輪和邊界質量塊開始進入同步旋轉，其展開過程的構形演化類似圖18。綜合分析得知，系統(tǒng)初始轉速對太陽帆展開至最大半徑所耗費的時間以及完全展開后發(fā)生的面外變形均會產生顯著影響。

圖25 不同初始轉速時邊界質量塊和中心彀輪的轉角Fig. 25 Rotation angle of edge mass and central hub at different initial speeds

圖26 不同初始轉速時的太陽帆展開半徑Fig. 26 Deployment radius of solar sail at different initial speeds

6 結 論

本文整合ANCF建模方法與張力場理論，針對薄膜結構拉壓不同模量的宏觀特性，提出了一種高效的建模與計算方法，為空間帆膜結構的展開動力學模擬提供了有效途徑。

1）針對薄膜太陽帆系統(tǒng)展開過程中顯著的非線性動力學特性，基于張力場理論與不同模量彈性理論，推導了單元精確的切向剛度矩陣，提出了一種考慮不同模量的ANCF薄膜單元計算框架。其中考慮了薄膜單元應力場對其本構關系的影響，準確地描述了薄膜結構在拉壓過渡區(qū)的非線性特性，有效提高了單元在牛頓迭代過程中的收斂性。

2）基于上述框架，采用第一類拉格朗日方程建立了六邊形薄膜太陽帆的自旋展開動力學模型，并通過可控數(shù)值耗散的廣義-α算法求解運動方程。根據數(shù)值模擬結果，分析了太陽帆不同結構設計參數(shù)對其展開book=228,ebook=132動力學特性的影響，為太陽帆類空間帆膜結構的系統(tǒng)設計和折展動力學仿真提供了一定的理論指導。

3）未來值得進一步研究的問題和方向包括：薄膜太陽帆完全展開之后，帆面產生的微弱面外變形以及振蕩過程仍未得到有效解決；現(xiàn)有的針對薄膜太陽帆自旋展開的動力學建模與仿真方法在計算效率和分析精度方面仍然存在不足，隨著薄膜單元網格的加密，難以針對大型帆膜結構進行高效模擬。