基于常應變單元的空間薄膜結構動力學建模與分析

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

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

2.西安電子科技大學電子裝備結構設計教育部重點實驗室,西安 710071

1 引言

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

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

本文由有限元離散方法將整個薄膜反射面離散為眾多三角形單元,這些單元的節點可視作質點,考慮這些三角形單元的彈性勢能,由最小勢能原理得到三角形單元的平衡方程,從而得到了薄膜反射面的有限元分析模型,進一步建立單元間的接觸模型,開展空間薄膜系統的展開動力學分析。基于薄板的下落算例,驗證了方法的正確性,隨后,在重力環境下,對薄膜反射面天線進行動力學展開分析,經過該展開動力學分析算例,對薄膜反射面展開過程的平穩性和動力學參數進行分析預測。

2 常應變三角形膜單元

空間薄膜結構往往構成拋物面反射面以滿足空間任務的需求,針對該類光滑連續的空間曲面,常基于有限單元法對目標曲面的變形予以建模。由有限元離散方法將整個薄膜反射面離散為眾多三角形單元,并建立對應拓撲關系的薄膜反射面數值分析模型,基于有限元建立了薄膜結構的展開動力學分析模型。為此,本節先建立用于離散空間薄膜結構的三角形單元,并給出該單元上點的應變和單元三條斜邊線應變之間的轉換關系。杜敬利等[21]提出三角形薄膜單元在不受力狀態時,可由其邊長描述膜單元的幾何構型。當三角形膜單元在不受力狀態下時,膜單元的幾何構型如圖1所示,該膜單元楊氏模量記作Em,薄膜厚度記作tm,泊松比記作vm。在三角形膜單元中,邊12、邊23及邊13所對應初始長度為l1、l2和l3。設定三角形膜單元的位移場為線性,由于應變為位移場的一階導數,故應變為常數,即為常應變三角形單元。當膜單元受力變形后,三角形膜單元的節點1、節點2和節點3對應坐標矢量則記為x1、x2和x3。膜單元受力拉伸后,沿邊12、邊23和邊13的應變如下:

(1)

應變的微分為:

(2)

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

ε=[εxεyγxy]T

(3)

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

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

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

(4)

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

所建立平面常應變三角形單元模型由三個節點坐標確定,單元上材料點的位置和變形是在全局坐標下定義的,簡化了幾何描述和坐標轉換,能夠滿足非線性收斂準則。

3 單元間的接觸模型

現任取薄膜反射面上的某單元123和對應膜單元的任意節點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為單元上任意點的坐標矢量。

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

(15)

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

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

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

(16)

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

(17)

且:

4 薄膜結構展開動力學方程

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

(18)

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

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

(19)

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

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

D=β1M+β2K

(20)

式中:β1和β2為比例系數。

5 算例分析

5.1 薄板單擺動力學仿真

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

圖3 薄板模型和網格分布

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

圖4 下落中的薄板構型和關鍵節點的位移曲線

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

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

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

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

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

圖6 剛性肋的運動學分析

(21)

節點n的速度矢量:

(22)

節點n的加速度矢量

(23)

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

圖7 展開角速度規劃

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

表1 反射面的物理參數

圖8 傘狀天線薄膜反射面網格模型

(24)

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

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

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

圖11 天線系統動能和應變能隨時間變化曲線

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

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

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

6 結論

本文構造了常應變三角形單元,建立了空間薄膜結構的動力學模型,通過薄板的單擺算例和傘狀薄膜反射面的展開過程算例,驗證了方法的準確性,并得到如下結論。

1)基于常應變三角形單元將薄膜反射面離散為數值曲面,根據離散的有限元單元和拓撲關系,能夠建立針對空間薄膜結構的動力學方程,能夠對大變形的薄膜結構動力學特性進行研究,由于膜本身具有的柔性特征,類似存在大變形膜結構的動力學展開模擬中,往往需要將步長設置的足夠小來保證數值收斂。

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

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