索網結構非線性有限元分析中初始張力狀態的快速構造方法

王 威，吳 松，郭其威，唐國安

（1.復旦大學 航空航天系，上海 200433；2.上海空間飛行器機構重點實驗室，上海 201108；3.上海宇航系統工程研究所，上海 201108）

0 引言

索網結構是張力結構體系的一個重要組成部分，其具有重量輕、強度高、收縮比大、成形多樣等顯著優點［1-2］，可被應用于太空可展開天線領域。預張力是索網結構成型的必要條件和提高其曲面穩定性的重要參數［3-4］。

在結構負荷工作階段，無論其靜力行為還是動力行為均體現出較高的非線性特性［5-6］，然而由于航天發展所需天線趨于大型化，結構全尺寸地面試驗難以開展［7］，所以在設計階段，有限元數值分析顯示出了巨大優勢。

有限元分析時，預張力的施加方法有等效荷載法、降溫法、初始應變法等［8-17］。而對于復雜的索網結構，等效載荷法并不適用，初應變法與降溫法本質相同，都是通過改變索網的單元原長度來施加預張力。根據是否進行迭代，確定初應變（或降溫量）的方法可分為迭代修正法和直接法。迭代修正法先取一個初應變（或降溫量）計算索單元的預張力，根據結果與設計值的差異修正初應變（或降溫量），然后再重新計算，直至迭代收斂［9-10，13］。在現階 段，這是一 種常用 的計算方法，但需要迭代計算，會導致效率降低。直接法不需要迭代，效率高。張航［11］和YANG 等［18］提出了考慮支座位移的計算方法，將總降溫值分為拉索兩端固定時達到設計預張力所需的降溫值，支撐結構位移帶來的拉索變形所對應的降溫值兩部分疊加。此種方法沒有考慮兩部分降溫值之間的耦合效應，應用于復雜的非線性結構時精度并不高。因此，需要一種能在利用有限元數值方法分析時，滿足高精度和高效率的初應變（或降溫量）的計算方法。

針對上述問題，本文綜合考慮索網的幾何非線性、剛度矩陣奇異性、框架的變形、預張力施加精度和計算效率等因素，提出了一種新的確定索網初應變（或降溫量）的方法。該方法以無預應力索網結構的有限元模型為基礎，迭加上帶預緊力柔索的初應力剛度矩陣，使得能夠借助通用程序MSC/NASTRAN 直接計算出降溫量，用于模擬索網的預張力。此方法綜合了迭代修正法和直接法的優勢，計算過程簡便，可操作性強，能保證計算精度，有效提高了效率。對應用于衛星可展開天線在內的復雜、多層、具有彈性約束的索網結構，沒有構型限制，且可進一步開展非線性固有頻率和響應的計算。

1 索單元預緊力降溫模擬的線性化計算方法

1.1 幾何非線性問題的單元切線剛度矩陣

由虛功原理可知，在平衡力系作用下的單元內，應力σ(e)在虛應變δε(e)上作的功等于等效節點力F(e)在虛位移δq(e)上作的功，即

根據連續介質力學的幾何方程，虛應變又可以用虛位移表示：

式中：BL為與變形無關的線性應變矩陣；BN為與位移q(e)有關的非線性應變矩陣。

將式（2）代入式（1）中，根據虛位移δq(e)的任意性，可以得到單元的平衡方程：

式中：Ω為體積域。

在材料屈服極限內，應力-應變關系依舊可以用胡克定律表示：

式中：D為材料的彈性矩陣；為單元的初應力向量；ε(e)為單元的應變向量。

將式（4）代入式（3）中，可得

因為矩陣BN和向量ε(e)均與q(e)有關，所以Ψ(e)是q(e)的非線性函數，式（5）是關于與q(e)的非線性方程組。非線性方程組一般采用迭代法求解，例如，具有二次收斂速度的牛頓-拉夫遜（Newton-Raphson，N-R）算法。N-R 算法需要提供函數Ψ(e)關于自變量q(e)的梯度

推導過程利用了關系?ε=BL+BN，這里符號?為梯度算子。若記

1.2 預張力索單元的初應力剛度矩陣

關于位移分量求變分后得到

對比式（10）與式（2）可知，索單元的線性和非線性應變矩陣為

再對式（12）的BN求變分，還能得到

式中：

式中：I3為3×3 的單位陣。

若給定索單元的預張力Ne，根據式（7）、式（8）以及式（11）、式（12）和式（14）便可計算出預張力索單元的初應力剛度矩陣

1.3 降溫值計算

對于由柔索和彈性構件組成的索網結構，用“節點-單元-材料-屬性”等常規的數據定義初始有限元模型。從模型數據中篩選出柔索單元，根據柔索張力的設計要求，按照式（15）逐個計算單元的初應力剛度矩陣Kσ(e)，再用常規的有限元裝配過程組裝成整個索網的初應力剛度矩陣Kσ。

將初應力剛度矩陣Kσ迭加到索網結構有限元模型的線性剛度矩陣上，近似地等效出索網結構帶有張緊力時的切線剛度矩陣。這樣處理后索網結構的剛度矩陣不再具有奇異性，而且是與變形或位移無關的常矩陣，可用線性方法進行計算。采用通用的有限元程序計算時，可以用降溫產生的收縮力等效柔索的預緊力。假設柔索單元的序號是1～nc，依次在柔索單元j上施加單位負溫度，計算出每一個柔索單元的張力，記為cij(i=1，2，…，nc)。當編號j從1 到nc遍歷全部柔索單元后，就得到了以cij為元素、溫度對張力的影響矩陣C，這時C是一個nc×nc非奇異階矩陣。如果柔索單元張力的設計值是N(i)(i=1，2，…，nc)，那么各個柔索單元上需要施加的負溫度T(i)(i=1，2，…，nc)應滿足關系

式中：N、T分別為由N(i)和T(i)排成的列向量。

上述方法在計算柔索單元的張力cij時，雖然略去了大變形剛度矩陣，但保留了切線剛度矩陣中占主導的初應力剛度矩陣和線性剛度矩陣，因此具有較高的計算精度。且在張力已知的情況下，初應力剛度矩陣是常矩陣，計算C矩陣的每一列可作為相同模型、不同工況問題求解，不需要重新生成和分解切線剛度矩陣，所以方法又具有較高的效率。

2 實施過程和簡例驗證

2.1 驗證模型的定義

用于說明算法的索網與結構組合體有限元模型的示意圖如圖1 所示。圖中，帶有#的數字為節點序號，小括號內的數字為單元序號。單元（1）～單元（5）是組成索網的五根柔索，材料為尼龍繩（彈性模量取1.5 GPA），其中，長 度l1=0.2 m 的橫索網（1）～橫索網（4）和長度l2=0.1 m 的豎索網（5）分別用截面抗拉剛度EA1=3 000 N 和EA2=150 N 的桿單元建模，支撐構件（6）和支撐構件（7）用剛度系數為k=105N/m 的接地彈簧建模。設計要求橫索網和豎索網的張力分別是N1=100 N，N2=1 N。考慮在平面內的變形，該索網結構組的 自由度為 6，對應的 位移分量u2、u3、u5、u6、v2、v5已標注圖1 中。

圖1 簡單模型的結構示意Fig.1 Structural sketch of simple model

算例用有限元分析程序NASTRAN 計算，模型的定義見表1。定義模型的數據包括節點坐標、單元組合信息、材料常數、單元屬性和約束條件。為了便于分析，將柔索材料的熱膨脹系數設置為α=1。采用這樣的設置，施加在柔索單元上溫度載荷T(e)就是單元單位長度的伸長量。由圖1可見，模型中節點2 和節點3 在縱向沒有剛度，索網與結構組合體的整體剛度矩陣是奇異的。為克服奇異性，在表1 中補充了第20 行Include‘KSIGMA.pch’，其內容是將要迭加的初應力剛度矩陣Kσ。

2.2 初應力剛度矩陣的生成

根據表1 中單元（1）～單元（5）的組合信息和相應的節點坐標，用算法語言自行編制程序，算得這5個單元構成的索網初應力剛度矩陣為

表1 索網結構算例的有限元模型定義Tab.1 Definitions of the finite element model for cable net structure

式中：Kσ的單位為N/m。

該矩陣是對稱的，其6 行（或者6 列）對應的位移分量在模型中分別是u2、v2、u3、u5、v5和u6。參照NASTRAN 的DMIG 格式要求，將Kσ命名為KSIGMA，下三角部分元素按行順序寫到外部文件，取名為KSIGMA.pch，見表2。

表2 初應力剛度矩陣的DMIG 數據Tab.2 DMIG data of the initial stress stiffness matrix

2.3 溫度對張力的影響矩陣生成

為了得到溫度對張力的影響矩陣C，再建立一個NASTRAN 的輸入文件，見表3。

表3 計算溫度對張力影響矩陣的NASTRAN 輸入文件Tab.3 NASTRAN input file for calculating the effect matrix of temperature on tension

該輸入文件包含3 個計算工況，分別對應于在柔索單元（1）、單元（2）和單元（5）上施加單位負溫度，求解5 個柔索單元的張力。由于影響矩陣以及結構均具有對稱性，單元（3）和單元（4）上施加單位負溫度的工況可以用工況1 和工況2 替代。運行NASTRAN 后，經整理得到

式中：矩陣中只有變量名、沒有具體數字的元素為可以利用結構或矩陣對稱性替代的數據。

2.4 耦合變形和張力的計算分析

根據五根柔索張力的設計值100、100、100、100和1 N，用式（17）給出的矩陣C根據式（16）可算得每個單元的溫度為

將這組溫度值作為已知單元溫度載荷條件，建立的NASTRAN 輸入文件，見表4。

表4 耦合變形和張力計算的NASTRAN 輸入文件Tab.4 NASTRAN input file for coupled deformation and tension calculation

計算后得到的柔索單元張力與設計值完全一致。關于變形的位移分量可以在微小位移假設下用解析法估算。參考圖1，節點3 的橫向位移就是在柔索（1）和柔索（2）張力N1作用下支撐彈簧（6）的伸長量，即柔索（1）和柔索（2）的收縮量相同，其和等于u3，因此-5×10-4m。節點2 豎向位移使得柔索（1）和柔索（2）的張力N1改變方向，其豎向的分量和與柔索（5）的張力N2平衡，由此得出-10-3m。位移分量的有限元計算與解析法估算結果的對比見表5，兩者相差無幾，表明用本文介紹的方法具有很高的計算精度。

表5 預緊力索網-結構算例耦合變形的計算結果Tab.5 Calculation results of the coupled deformation of the cable net structure with pretension

配置柔索等效溫度下降量、后續的耦合變形和張力計算都是線性過程，不需要迭代就能取得很高的計算效率。當預緊力的等效溫度確定后，作為已知載荷問題也可以用非線性分析計算索網-結構的耦合變形。此時，將橫索網和豎索網分別等效為邊的方型截面梁單元，利用細長梁較弱的橫向剛度避免零張力狀態下索網剛度矩陣的奇異性，實現非線性計算的初始迭代。用NASTRAN 求解序列106 計算，得到的橫索網張力均是100 N，與設計值完全吻合，豎索網張力是0.974 4 N，誤差約為-2.56%，符合工程計算的精度要求。位移的非線性計算結果見表5。解析法與線性和非線性方法得到的結果誤差相當。

3 索網式天線動力學特性計算的工程應用

索網式反射面天線如圖2 所示，結構包括伸展臂、支撐桁架、索網3 個主要部分。支撐桁架是由輔助驅動單元、外生支撐單元、立柱或助桿接頭組成，索網包含前索網、背索網、張緊索和斜拉索。天線支撐桁架左右跨度為14.76 m，前后跨度為12.00 m。

圖2 索網式可展開天線Fig.2 Cable net deployable antenna

運用上述算法，根據索網預張力設計值，配置索網單元的等效溫度下降量并作為輸入條件，然后采用SOL 106 求解序列，進行非線性靜力分析。分析時，設定伸展臂根部固支，索網預張力施加結果如圖3 所示，與設計值相比，相對誤差處于2.3%以下，因此能滿足設計階段的精度要求。在預張力作用下的索網、支撐桁架和伸展臂耦合變形的計算結果如圖4 所示。圖中，粗實線為變形放大50 倍后的狀態，節點位移的最大值在厘米量級。結合預張力施加精度及變形結果可知，以線性剛度矩陣和初應力剛度矩陣的疊加等效帶有張緊力的索網結構的切線剛度矩陣，在實際工程應用中是可行的。

圖3 索網預張力施加結果Fig.3 Result of pretension application to cable net

圖4 索網、框架和伸展臂的耦合變形Fig.4 Coupling deformation of cable net，frame，and extension arm

同樣采用SOL 106 求解序列，根據存在預張力狀態時的切線剛度矩陣進行模態分析，天線的前6階固有頻率以及與其對應的振型如圖5 所示。圖中：第1 階和第2 階主要為伸展臂的彎曲模態，由伸展臂自身剛度決定；第3 階為天線整體的1 階扭轉模態；第4 階主要為支撐桁架與索網結構的1 階扭轉模態；第5 階主要為伸展臂的局部模態，受到伸展臂自身連接剛度影響；第6 階為天線整體的2 階扭轉模態。

圖5 非線性固有模態Fig.5 Nonlinear natural mode

4 結束語

本文針對簡單柔索算例的詳細分析，以及在實際索網結構工程中的應用，表明對于幾何非線性效應顯著、橫向剛度中張緊力占主導的索網單元，在單元切線剛度矩陣中略去因位形改變而引起的初位移剛度矩陣，而保留由張力所產生的初應力剛度矩陣，能夠有效計算出用柔性索單元的預緊力。由于保留的初應力剛度矩陣是常量，因此，預緊力的計算屬于線性問題，具有很高的計算效率。將所得的預緊力用降低溫度的方法比擬，用熱彈性有限元計算出的柔索張力能夠與設計值高度吻合。對于航天器上應用的索網天線，本方法能在保證精度的前提下提高計算效率。此外，整個過程可以充分利用通用有限元計算程序，可以在此基礎上進一步開展模態和響應等力學計算，便于在工程中應用和實施。