基于常應變單元的空間薄膜結(jié)構(gòu)動力學建模與分析

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

1.湖南工學院 智能制造與機械工程學院,衡陽 421002

2.西安電子科技大學電子裝備結(jié)構(gòu)設計教育部重點實驗室,西安 710071

1 引言

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

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

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

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

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

(1)

應變的微分為:

(2)

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

ε=[εxεyγxy]T

(3)

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

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

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

(4)

由公式(4),進一步可聯(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é)點上的力向量可經(jīng)由外荷載P轉(zhuǎn)換得到,K(lt,xt)是關于節(jié)點位置矢量xt和單元邊長lt的總體剛度矩陣。

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

3 單元間的接觸模型

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

圖2 薄膜接觸判斷

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

(12)

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

(13)

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

進一步可得到平面的方程:

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

(14)

式中:x為單元上任意點的坐標矢量。

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

(15)

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

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

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

(16)

而當質(zhì)點A處于單元123平面的下方時:

(17)

且:

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

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

(18)

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

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

(19)

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

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

D=β1M+β2K

(20)

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

5 算例分析

5.1 薄板單擺動力學仿真

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

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

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

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

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

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

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

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

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

圖6 剛性肋的運動學分析

(21)

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

(22)

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

(23)

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

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

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

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

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

(24)

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

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

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

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

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

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

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

6 結(jié)論

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

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

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

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