杜雪林,張康,張瑞翔,葉拓,易文慧,蔣嘉斌

1.湖南工學(xué)院 智能制造與機(jī)械工程學(xué)院,衡陽 421002

2.西安電子科技大學(xué)電子裝備結(jié)構(gòu)設(shè)計教育部重點(diǎn)實(shí)驗(yàn)室,西安 710071

1 引言

目前,載人航天、深空探測、大型空間望遠(yuǎn)鏡和對地觀測等航天任務(wù)的開展和推進(jìn),無不對空間薄膜結(jié)構(gòu)的應(yīng)用提出了迫切要求。空間薄膜結(jié)構(gòu)質(zhì)量輕、易折疊,可被應(yīng)用在靜電成形薄膜反射面天線、充氣式可展開天線及太陽帆等航天器中,可作為新型空間展開結(jié)構(gòu)解決空間結(jié)構(gòu)日益增大及運(yùn)載火箭體積有限之間的矛盾[1-4]。為了提高星載天線的在軌性能和可靠性,需要進(jìn)一步系統(tǒng)開展采用大型可展開天線的動力學(xué)研究[5]。

空間膜結(jié)構(gòu)是由高強(qiáng)度膜材料通過一定方式施加預(yù)張力從而形成特定空間形狀并能承擔(dān)相應(yīng)外部載荷的一種結(jié)構(gòu)[6]。對薄膜結(jié)構(gòu)的動力學(xué)分析方面,往往采用大量的薄殼、索和膜等大變形和強(qiáng)幾何非線性單元,因此在展開過程中會發(fā)生大變形和大轉(zhuǎn)動,具有強(qiáng)非線性、復(fù)雜接觸等特點(diǎn),動力學(xué)建模與分析難度較大[7]。近年來,不少研究者對大型薄膜天線的動力學(xué)建模進(jìn)行了研究。如,Shen等[8]利用ABAQUS分析了張拉索提供的預(yù)張載荷發(fā)生變化時膜天線的模態(tài)振型和固有頻率[8]。Fan等研究了大型膜天線衛(wèi)星系統(tǒng)的動力學(xué)建模和在軌參數(shù)辨識問題,利用耦合系數(shù)法建立了衛(wèi)星系統(tǒng)的平動、轉(zhuǎn)動和振動動力學(xué)方程[9]。Yuan等研究了可折疊折紙空間膜結(jié)構(gòu)展開過程中考慮接觸沖擊的動力學(xué)建模問題[10]。余建新等建立了充氣展開反射面的有限元模型,采用氣囊模型和控制氣體法分析反射面從纏繞狀態(tài)到完全展開狀態(tài)的全過程[11]。王永等針對太陽帆航天器的薄膜帆面和大尺度桿的柔性特征,建立了相應(yīng)的非線性剛?cè)狁詈蟿恿W(xué)模型,并研究了對應(yīng)的控制方法[12]。榮吉利等采用有限元數(shù)值分析的方法,在 NASA 模態(tài)分析數(shù)值模型的基礎(chǔ)上,建立了圓形薄膜太陽翼展開動力學(xué)模型,以不同的展開控制函數(shù)為變量,研究了不同控制函數(shù)對展開平穩(wěn)性的影響[13]。Lampani等采用非線性有限元法對三維充氣結(jié)構(gòu)的展開機(jī)理進(jìn)行分析,指出當(dāng)問題的幾何變化不極端時,隱式有限元解可以提供更多的結(jié)果,而顯式方法能夠求解更復(fù)雜的充氣系統(tǒng)的時間解[14]。肖瀟等采用彈簧-質(zhì)點(diǎn)模型,對充氣管的充氣展開過程進(jìn)行數(shù)值模擬,得到在無重力環(huán)境下不同時刻卷曲折疊管的展開構(gòu)型[15]。Xu等對整個天線結(jié)構(gòu)進(jìn)行了參數(shù)化建模和分析,采用罰函數(shù)法求解薄膜的自接觸,基于改進(jìn)的彈簧-質(zhì)量系統(tǒng)對天線展開機(jī)構(gòu)進(jìn)行了仿真[16]。Borse等通過在MATLAB中進(jìn)行有限元應(yīng)力分析,對方形、扁平、薄充氣膜結(jié)構(gòu)的性能進(jìn)行預(yù)測[17]。Li等建立了包含索單元、膜單元和梁單元的集成有限元模型,對鋼支撐結(jié)構(gòu)的張拉膜進(jìn)行分析[18]。Liu等建立了具有參數(shù)不確定性的阻尼膜結(jié)構(gòu)動力分析的隨機(jī)動力剛度解析公式,所開發(fā)的隨機(jī)動力剛度單元可以在考慮不確定性的一般邊界條件下組裝成膜組件的模型[19]。Liu等提出了一種將褶皺/松弛模型集成到柔性多體動力學(xué)絕對節(jié)點(diǎn)坐標(biāo)系薄殼單元中的計算方法,研究展開太陽帆和充氣氣囊等膜系統(tǒng)動力學(xué)分析中的褶皺現(xiàn)象[20]?？偟膩砜?目前空間薄膜結(jié)構(gòu)的研究主要采用增量有限元方法,由于膜結(jié)構(gòu)的大位移和大變形,接觸等因素的,往往導(dǎo)致剛度矩陣奇異,造成求解困難。

本文由有限元離散方法將整個薄膜反射面離散為眾多三角形單元,這些單元的節(jié)點(diǎn)可視作質(zhì)點(diǎn),考慮這些三角形單元的彈性勢能,由最小勢能原理得到三角形單元的平衡方程,從而得到了薄膜反射面的有限元分析模型,進(jìn)一步建立單元間的接觸模型,開展空間薄膜系統(tǒng)的展開動力學(xué)分析?；诒“宓南侣渌憷?驗(yàn)證了方法的正確性,隨后,在重力環(huán)境下,對薄膜反射面天線進(jìn)行動力學(xué)展開分析,經(jīng)過該展開動力學(xué)分析算例,對薄膜反射面展開過程的平穩(wěn)性和動力學(xué)參數(shù)進(jìn)行分析預(yù)測。

2 常應(yīng)變?nèi)切文卧?/h2>

空間薄膜結(jié)構(gòu)往往構(gòu)成拋物面反射面以滿足空間任務(wù)的需求,針對該類光滑連續(xù)的空間曲面,?；谟邢迒卧▽δ繕?biāo)曲面的變形予以建模。由有限元離散方法將整個薄膜反射面離散為眾多三角形單元,并建立對應(yīng)拓?fù)潢P(guān)系的薄膜反射面數(shù)值分析模型,基于有限元建立了薄膜結(jié)構(gòu)的展開動力學(xué)分析模型。為此,本節(jié)先建立用于離散空間薄膜結(jié)構(gòu)的三角形單元,并給出該單元上點(diǎn)的應(yīng)變和單元三條斜邊線應(yīng)變之間的轉(zhuǎn)換關(guān)系。杜敬利等[21]提出三角形薄膜單元在不受力狀態(tài)時,可由其邊長描述膜單元的幾何構(gòu)型。當(dāng)三角形膜單元在不受力狀態(tài)下時,膜單元的幾何構(gòu)型如圖1所示,該膜單元楊氏模量記作Em,薄膜厚度記作tm,泊松比記作vm。在三角形膜單元中,邊12、邊23及邊13所對應(yīng)初始長度為l1、l2和l3。設(shè)定三角形膜單元的位移場為線性,由于應(yīng)變?yōu)槲灰茍龅囊浑A導(dǎo)數(shù),故應(yīng)變?yōu)槌?shù),即為常應(yīng)變?nèi)切螁卧?。?dāng)膜單元受力變形后,三角形膜單元的節(jié)點(diǎn)1、節(jié)點(diǎn)2和節(jié)點(diǎn)3對應(yīng)坐標(biāo)矢量則記為x1、x2和x3。膜單元受力拉伸后,沿邊12、邊23和邊13的應(yīng)變?nèi)缦?

(1)

應(yīng)變的微分為:

(2)

在圖1中,選擇三角形膜單元節(jié)點(diǎn)1為坐標(biāo)系原點(diǎn),將邊12方向看作x軸,取坐標(biāo)系平面和平面123重合,即得到局部坐標(biāo)系Omxmym。在局部坐標(biāo)系中,任意位于該三角形上的點(diǎn)的應(yīng)變表示:

ε=[εxεyγxy]T

(3)

式中:εx和εy沿xm和ym方向的正應(yīng)變記作εx和εy,切應(yīng)變記為γxy。

現(xiàn)將xm軸由當(dāng)前位置逆時針轉(zhuǎn)動,旋轉(zhuǎn)角度為α0,沿該方向的正應(yīng)變記為εp,基于公式(3),可得到εp的表達(dá)式:

εp=cos2α0εx+sin2α0εy+cosα0sinα0γxy

(4)

由公式(4),進(jìn)一步可聯(lián)系公式(1)和(3),得到:

εl=Ttε

(5)

式中:

三角形膜單元上的彈性勢能和外力做功之和為:

(6)

Am=[p(p-l1)(p-l2)(p-l3)]1/2

(7)

基于式(1)、(2)和(5),勢能的微分為:

(8)

上式可另寫作:

(9)

因此,基于最小勢能原理,得到了三角形單元的平衡方程為:

km(lm,xm)xm-fm=0

(10)

K(lt,xt)xt=F

(11)

式中:F為節(jié)點(diǎn)上的力向量可經(jīng)由外荷載P轉(zhuǎn)換得到,K(lt,xt)是關(guān)于節(jié)點(diǎn)位置矢量xt和單元邊長lt的總體剛度矩陣。

所建立平面常應(yīng)變?nèi)切螁卧Ｐ陀扇齻€節(jié)點(diǎn)坐標(biāo)確定,單元上材料點(diǎn)的位置和變形是在全局坐標(biāo)下定義的,簡化了幾何描述和坐標(biāo)轉(zhuǎn)換,能夠滿足非線性收斂準(zhǔn)則。

3 單元間的接觸模型

現(xiàn)任取薄膜反射面上的某單元123和對應(yīng)膜單元的任意節(jié)點(diǎn)A,如圖2所示,可知單元123上的邊21、31和32的方向矢量寫為:

圖2 薄膜接觸判斷

n21=x2-x1,n31=x3-x1,n32=x3-x2

(12)

該薄膜單元所在平面的法向矢量可記作:

(13)

其中nV為該平面的法向單位矢量。

進(jìn)一步可得到平面的方程:

nV(x-xi)=0,i=?(1,2,3)

(14)

式中:x為單元上任意點(diǎn)的坐標(biāo)矢量。

現(xiàn)需要判斷點(diǎn)A相對于平面的位置,有nV·(xA-xi)>0,i=?(1,2,3)時質(zhì)點(diǎn)A位于平面123的上方,否則,則處于平面123的下方,如圖2所示,質(zhì)點(diǎn)A此時位于平面的上方。

(15)

可通過對投影點(diǎn)A′和三角形節(jié)點(diǎn)1、2和3分別形成的三角形A′12、A′13和A′23面積和來判別A′是否位于三角形單元上。即,如SΔA′12+SΔA′13+SΔA′23=SΔ123時,投影點(diǎn)處在三角形單元區(qū)域內(nèi),反之,投影點(diǎn)A′處于三角形單元區(qū)域外。

因此,可知如果投影點(diǎn)A′位于三角形單元123上,則點(diǎn)A和三角形單元123構(gòu)成接觸對,存在接觸的可能性。此時,需要確保點(diǎn)A和單元123之間的距離不得小于給定的距離,否則就需要施加懲罰力,避免接觸的發(fā)生。

在此,給定懲罰力Fp,如果點(diǎn)A位于三角形單元123平面的上方時:

(16)

而當(dāng)質(zhì)點(diǎn)A處于單元123平面的下方時:

(17)

且:

4 薄膜結(jié)構(gòu)展開動力學(xué)方程

空間薄膜結(jié)構(gòu)可被離散為眾多的三角形單元,可將這些三角形單元的節(jié)點(diǎn)看作對應(yīng)質(zhì)點(diǎn),對應(yīng)質(zhì)點(diǎn)上的質(zhì)量與其所處于的三角形的形狀和數(shù)目相關(guān),由三角形的面積求質(zhì)點(diǎn)的質(zhì)量,即:

(18)

式中:ρM為膜材料的密度;Ak為連接質(zhì)點(diǎn)i的第k個三角形單元的面積;s為連接質(zhì)點(diǎn)i的所有三角形單元數(shù)量;αi,k為以質(zhì)點(diǎn)i做頂點(diǎn)的第k個三角形單元的內(nèi)角。

從而可知,整個空間薄膜反射面可經(jīng)離散為眾多的三角形單元,把得到的節(jié)點(diǎn)看作質(zhì)點(diǎn),根據(jù)對應(yīng)的單元間拓?fù)潢P(guān)系,組集整個薄膜反射面的質(zhì)量矩陣和剛度矩陣,得到最終的展開動力學(xué)方程:

(19)

式中:M、D和K分別為薄膜反射面的總體質(zhì)量矩陣、比例阻尼矩陣和總體剛度矩陣;FI為節(jié)點(diǎn)上的慣性力;FE為系統(tǒng)所受外力矢量;Fp為系統(tǒng)界面需引入的懲罰力。

采用Rayleigh比例阻尼矩陣[22],則D可寫為:

D=β1M+β2K

(20)

式中:β1和β2為比例系數(shù)。

5 算例分析

5.1 薄板單擺動力學(xué)仿真

由動力學(xué)模型式(20),先建立重力環(huán)境下的空間薄板單擺模型,進(jìn)行算例分析,此處薄板參數(shù)同文獻(xiàn)[23],空間薄板的幾何參數(shù)如圖3(a)所示,板材料彈性模量取0.1MPa,材料密度為7810kg/m3,板材厚度為0.01m,泊松比取0.3。在薄板A點(diǎn)有一固定球鉸,在重力的作用下薄板做單擺運(yùn)動。采用三角形單元對薄板進(jìn)行網(wǎng)格劃分,并標(biāo)示了一些關(guān)鍵節(jié)點(diǎn),在圖3(b)中,網(wǎng)格模型中,共劃分三角形單元200個,節(jié)點(diǎn)數(shù)目121。設(shè)定單擺時間為0.3s,進(jìn)行動力學(xué)仿真分析。

圖3 薄板模型和網(wǎng)格分布

經(jīng)過計算得到如下結(jié)果,不同時間對應(yīng)的薄板構(gòu)型如圖4(a)所示,結(jié)果同文獻(xiàn)[23]基本一致。各關(guān)鍵點(diǎn)的位移曲線,即這些節(jié)點(diǎn)沿x、y和z軸的位移曲線分別示于圖4(b-d)。

圖4 下落中的薄板構(gòu)型和關(guān)鍵節(jié)點(diǎn)的位移曲線

由圖可見,處于同邊的節(jié)點(diǎn)的位移在下落前期運(yùn)動狀態(tài)一致,如P4和P7、P2和P3等,位置對稱的節(jié)點(diǎn)如P4和P2、P7和P3、P6和P8等其位移變化也具有對稱性(如P4、P7和P8沿x方向的位移曲線分別同P2、P3和P6沿y方向的位移曲線一致,反之亦然)。從結(jié)果上看,本仿真分析算例得到的結(jié)果合理且準(zhǔn)確。

5.2 徑向肋薄膜傘狀天線的展開動力學(xué)分析

徑向肋薄膜傘狀天線具有較為簡單的結(jié)構(gòu)形式,但是可靠性強(qiáng)并且展開后表面精度高,展開機(jī)制和雨傘類似,因此應(yīng)用廣泛。徑向肋薄膜傘狀天線主要有徑向肋和薄膜反射面構(gòu)成,在徑向肋之間連接有剪裁的扇形薄膜反射面,從而形成整個拋物反射面。傘狀天線的收展運(yùn)動是由徑向肋的運(yùn)動帶動薄膜反射面實(shí)現(xiàn)的,在此,給出了共6個剛性肋的徑向肋薄膜傘狀天線的簡化模型,如圖5所示,隨后對天線的展開過程進(jìn)行動力學(xué)分析。

圖5 簡化的傘狀天線反射面

將高強(qiáng)度的徑向肋看作剛體,先結(jié)合徑向肋的幾何性質(zhì)作運(yùn)動學(xué)分析,任意取一根剛性肋,如圖6所示。將選定的剛性肋離散為n-1個單元,共n個節(jié)點(diǎn),在剛性肋的收展運(yùn)動中,各個節(jié)點(diǎn)的位移可經(jīng)由幾何關(guān)系進(jìn)行推導(dǎo)。

圖6 剛性肋的運(yùn)動學(xué)分析

(21)

節(jié)點(diǎn)n的速度矢量:

(22)

節(jié)點(diǎn)n的加速度矢量

(23)

對剛性肋進(jìn)行運(yùn)動學(xué)分析,設(shè)定圖5中的天線孔徑為0.55m,拋物面焦距1.1m,展開時間設(shè)定為1s,展開角度α=60°,對角速度進(jìn)行規(guī)劃,規(guī)劃曲線見圖7。在t=0s時天線處于收攏狀態(tài),而t=1s時天線完全展開。

圖7 展開角速度規(guī)劃

綜合考慮剛性肋和薄膜反射面,建立徑向肋薄膜傘狀天線整體動力學(xué)模型。薄膜反射面的物理參數(shù),如表1所示,膜材料的密度為1390kg/m3,膜的現(xiàn)時構(gòu)型內(nèi)預(yù)應(yīng)變狀態(tài)為膜內(nèi)各點(diǎn)沿徑向應(yīng)變1.338×10-4,緯向應(yīng)變1.338×10-4,阻尼系數(shù)β1和β2分別取0.5 和0.001。剛性肋的運(yùn)動帶動薄膜反射面收展,由于柔性薄膜構(gòu)型未知,在此,先由完全展開態(tài)構(gòu)型逆向收攏得到傘狀天線的收攏態(tài)構(gòu)型,再逐漸由收攏態(tài)進(jìn)行展開分析。因此,把整個薄膜反射面離散為1200個三角形單元,共計661個節(jié)點(diǎn),網(wǎng)格劃分模型圖如圖8(a),該模型為軸對稱模型,因此可取1/6模型上三個關(guān)鍵節(jié)點(diǎn)跟蹤位移和速度,關(guān)鍵節(jié)點(diǎn)的位置和編號如圖8(b)。在展開過程中,由于反射面的中心和剛性肋中心鉸處固連,為此,需要施加對應(yīng)的位移約束:

表1 反射面的物理參數(shù)

圖8 傘狀天線薄膜反射面網(wǎng)格模型

(24)

展開過程中傘狀天線的構(gòu)型圖如圖9,由圖中可以觀察到展開過程中反射面的構(gòu)型變化。取展開時間為t=0.1s、t=0.3s、t=0.7s和t=0.9s,提取對應(yīng)時刻的薄膜反射面的位移,繪制位移云圖,如圖10所示,進(jìn)一步可明確展開過程中薄膜反射面的變形情況。由圖可知,整個展開過程,薄膜上節(jié)點(diǎn)都表現(xiàn)出良好的對稱性,節(jié)點(diǎn)位移由外徑向內(nèi)徑呈逐漸減小趨勢,最大節(jié)點(diǎn)位移發(fā)生在剛性肋頂點(diǎn)處,同實(shí)際情況一致。

圖10 薄膜反射面的位移云圖

為了清晰表征展開性能,求解并繪制了整個過程的薄膜反射面動能和應(yīng)變能變化曲線,如圖11所示,由圖11(a)觀之,0～0.4s 區(qū)間反射面上節(jié)點(diǎn)速度變化幅度大,然后漸趨平穩(wěn),后期剛性肋展至設(shè)計位置后,動能漸趨于零?？傮w來看,動能曲線存在局部震蕩,也說明了膜本身具有的柔性特征。需要指出,類似存在大變形膜結(jié)構(gòu)的動力學(xué)展開模擬中,往往需要將步長設(shè)置的足夠小來保證數(shù)值收斂。由圖11(b)觀之,約在0～0.8s區(qū)間,應(yīng)變能曲線整體平穩(wěn)但存在振蕩,這期間肋間薄膜未被完全張緊,在0.8-1s區(qū)間,整個薄膜反射面在剛性肋的張緊作用下,逐漸完全展開,對應(yīng)地觀察到應(yīng)變能曲線在該區(qū)間呈上升趨勢。

圖11 天線系統(tǒng)動能和應(yīng)變能隨時間變化曲線

薄膜反射面上關(guān)鍵節(jié)點(diǎn)的位移曲線能夠表征系統(tǒng)展開過程的動力學(xué)響應(yīng),這些節(jié)點(diǎn)的位移變化情況見圖12(a～c)所示,如從位移曲線來看,節(jié)點(diǎn)P1、P2和P3位移曲線走勢平穩(wěn),可知展開過程較為順利且平穩(wěn)性好。觀察得知,P1沿三個方向的位移均大于其他兩個節(jié)點(diǎn)位移,這是因?yàn)镻1位于反射面最外環(huán)和肋的中間位置,變形最大。進(jìn)一步地,分析關(guān)鍵節(jié)點(diǎn)的速度變化趨勢,速度變化情況見圖12(d～f),觀察發(fā)現(xiàn),速度曲線存在明顯的震蕩現(xiàn)象,由幅度來看,圖12(d)中,P3處速度變化幅度較大,

圖12 天線關(guān)鍵節(jié)點(diǎn)位移和速度隨時間變化曲線

這是由于該點(diǎn)位移曲線變化幅度大所導(dǎo)致的。以此類推,在圖12(e～f)中,P1點(diǎn)處速度變化幅度較大。圖12(e～f)中節(jié)點(diǎn)沿z軸的速度曲線的振動幅度最大,由于膜結(jié)構(gòu)自身的顯著柔性,可知薄膜的位移變形主要是沿z軸方向發(fā)生,因此對z軸方向進(jìn)行振動控制尤為必要。

6 結(jié)論

本文構(gòu)造了常應(yīng)變?nèi)切螁卧?建立了空間薄膜結(jié)構(gòu)的動力學(xué)模型,通過薄板的單擺算例和傘狀薄膜反射面的展開過程算例,驗(yàn)證了方法的準(zhǔn)確性,并得到如下結(jié)論。

1)基于常應(yīng)變?nèi)切螁卧獙⒈∧し瓷涿骐x散為數(shù)值曲面,根據(jù)離散的有限元單元和拓?fù)潢P(guān)系,能夠建立針對空間薄膜結(jié)構(gòu)的動力學(xué)方程,能夠?qū)Υ笞冃蔚谋∧そY(jié)構(gòu)動力學(xué)特性進(jìn)行研究,由于膜本身具有的柔性特征,類似存在大變形膜結(jié)構(gòu)的動力學(xué)展開模擬中,往往需要將步長設(shè)置的足夠小來保證數(shù)值收斂。

2)通過求解并繪制了整個過程的薄膜反射面動能和應(yīng)變能變化曲線,能夠表征系統(tǒng)展開過程的動力學(xué)響應(yīng),分析展開的平穩(wěn)性。需要特別指出的是,本文方法還可以靈活擴(kuò)展,通過替換膜單元彈性系數(shù)矩陣,可實(shí)現(xiàn)采用正交異性膜材的膜反射面進(jìn)行展開動力學(xué)分析。

實(shí)際上,薄膜結(jié)構(gòu)的展開過程十分復(fù)雜,如何將更多的影響因素考慮進(jìn)來,將展開模型建立的更為準(zhǔn)確,如更為貼合空間薄膜結(jié)構(gòu)的折紙技術(shù)的應(yīng)用,討論膜結(jié)構(gòu)褶皺和松弛的影響也至關(guān)重要,特別是褶皺松弛造成的薄膜剛度矩陣的奇異使求解變得更加困難。基于該思路,進(jìn)一步優(yōu)化展開動力學(xué)模型,是下一步研究的問題。