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

（1.復(fù)旦大學(xué) 航空航天系，上海 200433；2.上海空間飛行器機(jī)構(gòu)重點(diǎn)實(shí)驗(yàn)室，上海 201108；3.上海宇航系統(tǒng)工程研究所，上海 201108）

0 引言

索網(wǎng)結(jié)構(gòu)是張力結(jié)構(gòu)體系的一個重要組成部分，其具有重量輕、強(qiáng)度高、收縮比大、成形多樣等顯著優(yōu)點(diǎn)［1-2］，可被應(yīng)用于太空可展開天線領(lǐng)域。預(yù)張力是索網(wǎng)結(jié)構(gòu)成型的必要條件和提高其曲面穩(wěn)定性的重要參數(shù)［3-4］。

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

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

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

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

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

由虛功原理可知，在平衡力系作用下的單元內(nèi)，應(yīng)力σ(e)在虛應(yīng)變δε(e)上作的功等于等效節(jié)點(diǎn)力F(e)在虛位移δq(e)上作的功，即

根據(jù)連續(xù)介質(zhì)力學(xué)的幾何方程，虛應(yīng)變又可以用虛位移表示：

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

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

式中：Ω為體積域。

在材料屈服極限內(nèi)，應(yīng)力-應(yīng)變關(guān)系依舊可以用胡克定律表示：

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

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

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

推導(dǎo)過程利用了關(guān)系?ε=BL+BN，這里符號?為梯度算子。若記

1.2 預(yù)張力索單元的初應(yīng)力剛度矩陣

關(guān)于位移分量求變分后得到

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

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

式中：

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

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

1.3 降溫值計算

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

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

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

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

2 實(shí)施過程和簡例驗(yàn)證

2.1 驗(yàn)證模型的定義

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

圖1 簡單模型的結(jié)構(gòu)示意Fig.1 Structural sketch of simple model

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

2.2 初應(yīng)力剛度矩陣的生成

根據(jù)表1 中單元（1）～單元（5）的組合信息和相應(yīng)的節(jié)點(diǎn)坐標(biāo)，用算法語言自行編制程序，算得這5個單元構(gòu)成的索網(wǎng)初應(yīng)力剛度矩陣為

表1 索網(wǎng)結(jié)構(gòu)算例的有限元模型定義Tab.1 Definitions of the finite element model for cable net structure

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

該矩陣是對稱的，其6 行（或者6 列）對應(yīng)的位移分量在模型中分別是u2、v2、u3、u5、v5和u6。參照NASTRAN 的DMIG 格式要求，將Kσ命名為KSIGMA，下三角部分元素按行順序?qū)懙酵獠课募∶麨镵SIGMA.pch，見表2。

表2 初應(yīng)力剛度矩陣的DMIG 數(shù)據(jù)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 個計算工況，分別對應(yīng)于在柔索單元（1）、單元（2）和單元（5）上施加單位負(fù)溫度，求解5 個柔索單元的張力。由于影響矩陣以及結(jié)構(gòu)均具有對稱性，單元（3）和單元（4）上施加單位負(fù)溫度的工況可以用工況1 和工況2 替代。運(yùn)行NASTRAN 后，經(jīng)整理得到

式中：矩陣中只有變量名、沒有具體數(shù)字的元素為可以利用結(jié)構(gòu)或矩陣對稱性替代的數(shù)據(jù)。

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

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

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

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

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

表5 預(yù)緊力索網(wǎng)-結(jié)構(gòu)算例耦合變形的計算結(jié)果Tab.5 Calculation results of the coupled deformation of the cable net structure with pretension

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

3 索網(wǎng)式天線動力學(xué)特性計算的工程應(yīng)用

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

圖2 索網(wǎng)式可展開天線Fig.2 Cable net deployable antenna

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

圖3 索網(wǎng)預(yù)張力施加結(jié)果Fig.3 Result of pretension application to cable net

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

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

圖5 非線性固有模態(tài)Fig.5 Nonlinear natural mode

4 結(jié)束語

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