大型空間結構熱-動力學耦合分析方法綜述1)

薛明德 向志海

(清華大學航天航空學院,北京 100084)

引言

航天器中有很多可展開附件,比如太陽翼、天線、桅桿等等,這些結構大多由薄壁桿件構成.隨著我國航天事業的進步,這類結構逐漸向大型化和柔性化方向發展.在軌展開之后,它們在太空失重環境中就主要受到熱載荷的作用,可能產生影響結構功能的熱致響應.比如準靜態熱變形會改變天線的形狀,導致其聚焦能力下降.而本文主要關注發生在低軌航天器上的熱誘發振動(thermally induced vibration)現象.因為低軌的地球半影區很窄,所以航天器只需要幾秒鐘就能進出地球陰影,從而使這類輕柔結構受到急劇變化的太陽熱流作用.對于這類低熱容和低剛度的結構,快速變化的熱流會迅速改變其截面溫差,并進一步激發振動響應,甚至產生不穩定的熱誘發振動,即熱顫振(thermal flutter)現象.

熱誘發振動現象在人類早期的航天活動中就已經出現.美國的第一顆衛星Explorer I 于1958 年發射入軌后很快就發生了姿態偏轉,其原因就是柔性天線的熱誘發振動耗散了能量[1].加拿大于1962 年發射的第一顆衛星Alouette 1 也出現了類似的事故[2].隨著航天事業的發展,熱誘發振動現象也隨著各種大型柔性結構的使用而陸續被觀測到[3],甚至出現過熱顫振現象.其中最著名的是1990 年發射的哈勃太空望遠鏡的柔性太陽翼彎扭耦合的熱顫振事故[4],直接影響了望遠鏡的成像精度.后來人們只能借助航天飛機將其太陽翼更換成半剛性結構,才使望遠鏡能正常工作.經過多次在軌維護后,目前的哈勃太空望遠鏡使用的是剛性太陽翼[5].哈勃熱顫振事故對國際航天工程產生了深遠的影響[6],對我國相關航天器的研制也具有重要的借鑒價值.

其實在人類開展航天活動之前,Boley[7]已從理論上預測了熱誘發振動現象,并引入了一個無量綱參數B來衡量發生彎曲型熱誘發振動時的動態位移dmax和穩態位移dst之間的比值,即熱誘發振動的劇烈程度[8]

式中,ω是結構的第一階固有振動圓頻率,代表了機械振動的快慢程度;τ稱為結構的熱特征時間,代表了截面溫差隨時間變化的快慢程度.該公式簡潔明了地解釋了為什么輕柔結構容易出現熱誘發振動.參數B也因此被稱為“Boley 參數”[9].不過式(1)只適用于突加熱流所產生的桿件彎曲型熱誘發振動,本文的第3 節將給出復雜的空間薄壁桿系結構任意形式熱誘發振動劇烈程度的度量公式,并討論它在漸變熱流作用下的準確性.

熱顫振機理的研究最早見于Augusti[10]對于一根同時受到軸力和熱流作用的桿模型的分析工作.他指出熱顫振是一種熱-對結構[11]強耦合效應所產生的自激振動現象.Donohue 和Frisch[11]也通過討論開口薄壁桿件的扭轉振動解釋了OV1-10 衛星桅桿的扭轉熱顫振現象.Yu[12]曾試圖建立懸臂桿彎曲熱顫振的判斷準則: 當變形前的懸臂桿指向太陽時就會發生熱顫振.但Graham[13]發現Yu[12]論文中的邊界條件有誤,并將該熱顫振準則修正為: 當變形前的懸臂桿背離太陽時就會發生熱顫振.該準則也被Thornton 和Kim[4]用來解釋哈勃太空望遠鏡的熱顫振事故.不過,人們在實驗中發現當熱流垂直于懸臂桿入射時也會發生熱顫振[3,14],并不符合Graham準則.造成理論預測和實驗結果不一致的原因是: 穩定性是非線性問題,必須針對變形后的平衡狀態進行分析,而Graham 熱顫振準則卻是基于變形前構型建立的.于是文獻[15]將Graham 熱顫振準則修正為: 當變形后的懸臂桿背離太陽時就會發生熱顫振.本文的第4 節將對此準則進行具體的介紹,并進一步推廣至復雜工程結構的熱顫振分析.

明確了物理機理后,人們更關心如何分析復雜工程結構的熱誘發振動和熱顫振響應.由于實際航天結構大多由開口和閉口薄壁桿件構成復雜的空間桿系,需要一體地分析該復雜結構中的傳熱學與動力學及其互相耦合問題,這種分析工作往往需要借助于有限元程序才能開展.Mason[16]于1968 年就開始用有限元程序計算簡支桿和板的熱誘發振動響應.但他的有限元模型中直接采用了溫度場的解析解,回避了輻射換熱條件下的瞬態溫度場分析這個高度非線性的難題.由于考慮輻射換熱的溫度計算工作量很大,用當時的計算機很難分析復雜結構的熱誘發振動問題.該狀況一直持續到20 世紀90 年代,特別是哈勃太空望遠鏡的熱顫振事故之后,人們才更加重視相關的研究.Namburu 和Tamma[17]直接采用三維有限元進行熱誘發振動分析,導致計算工作量非常大.為了提高計算效率,為文獻[18]將閉口薄壁桿的橫截面溫度場分解為平均溫度和攝動溫度兩部分,從而發展出一種譜單元來分析桿件橫截面的穩態溫度場,并可求出結構的準靜態熱變形.不過這種單元的平均溫度和攝動溫度相互耦合,計算工作量也不小.Thornton 和Kim[4]采用Ritz 法分析哈勃太空望遠鏡太陽翼的熱誘發振動問題時,將薄壁圓管桿件橫截面的溫度場表示成平均溫度加余弦攝動溫度的形式,簡化了熱傳導控制方程.而文獻[19-20]進一步將任意形式(閉口或單支開口)薄壁桿件橫截面的溫度場展開成一般的Fourier 級數形式,發展出一種全新的一維薄壁桿件單元.該單元的每一個節點包含一個平均溫度自由度和若干攝動溫度自由度.利用Fourier 級數的正交性,可將瞬態熱傳導方程解耦為僅包含平均溫度的非線性方程組和同時包含平均溫度和攝動溫度的線性方程組.因為有限元模型中的平均溫度自由度少,因此可以很快地求解該非線性方程組,然后再高效地求解含有攝動溫度自由度的線性方程組,從而提高了計算效率.另外,采用該瞬態溫度單元可以與結構動力學分析的桿件單元共用節點,因此非常有利于進行熱-結構的耦合分析[20-21].本文的第1 節將詳細介紹Fourier 溫度單元,并分析了平均溫度和攝動溫度時變特性的差別.本文的第2 節將介紹基于Fourier 溫度單元的熱誘發振動分析方法,主要強調在開口薄壁桿件、復雜結構中瞬態熱傳導-動力學耦合和幾何非線性分析方面需要注意的問題.

到20 世紀末,國際上對熱誘發振動的研究已經比較完整[22],但是分析精度尚不能完全滿足蓬勃發展的太空望遠鏡等具有高指向精度要求的航天器的需要[6,23].近年來,NASA 還進一步開展了在軌熱誘發振動實驗[24],相關成果已經應用于2021 年的“雙小行星重定向測試 (Double Asteroid Redirection Test,DART)”任務[25].我國的相關研究始于1997 年的國家高技術研究發展計劃.和其他國家的研究相比,雖然起步很晚,但是一直立足國情致力于發展能解決復雜工程問題的高效分析方法,并在熱顫振準則方面做出了有特色的工作.該研究工作斷斷續續地持續了二十多年,最近幾年才因我國航天工程的迫切需求而引起重視,并通過地面測試驗證了所發展的理論和方法[26-29],也在工程中得到了應用.在本文的后續小節中,將詳細地總結這個領域的研究成果,并展望未來的工作,以期為我國航天事業的發展做一些微薄的貢獻.限于篇幅,本文只強調解決這種強非線性問題(輻射換熱、幾何大變形、熱-動力學耦合)的思路,而略去了很多具體的數學推導.有興趣的讀者請在參考文獻中尋找更多的細節.

1 輻射換熱下薄壁桿件的瞬態溫度場分析

航天器可展附件通常包含由許多薄壁桿件組成的復雜構件.它們在太空中受到太陽和地球熱流照射,同時向空間輻射散熱.研究這類結構的熱誘發振動問題,就必須發展高效的瞬態溫度場分析算法,以及能與之協調工作的結構動力學分析算法.雖然有限元方法特別適合分析這種復雜結構,但是傳統溫度桿件單元的節點自由度只有平均溫度,無法刻畫桿件橫截面內的溫度分布,從而無法算出對熱誘發振動至關重要的熱彎矩和熱雙力矩.采用溫度殼單元或溫度實體單元雖然可以解決上述問題[17,30-31],但是需要在薄壁桿件橫截面內劃分很多單元,從而在求解輻射換熱這種強非線性方程時會花費較多的計算工作量.另外,這些溫度單元和結構單元往往不能共用一套有限元網格,不利于分析復雜工程結構的熱-動力學耦合響應.針對這些問題,提出了一種Fourier 薄壁桿件溫度有限單元,通過增加節點攝動溫度的方法來刻畫桿件橫截面的溫度分布,并利用Fourier 級數的正交性來提高計算效率.該溫度單元可以和結構動力學分析的桿單元共用節點,因此采用一套有限元網格就可以高效而準確地計算出復雜工程結構的熱誘發振動響應.

1.1 薄壁桿件溫度單元[19-20]

如圖1 所示,薄壁桿件單元的幾何尺寸中包含軸線長度l、橫截面周長Sh和壁厚ts三個尺度.由于壁厚遠小于其他兩個尺度,所以可以忽略溫度沿壁厚的變化.于是在分析隨時間t變化的溫度場T時只考慮它在桿件軸線和截面周向的變化,即T=T(θ,ξ,t),其中ξ=x/l(0 ≤x≤l) 和θ=2πs/Sh(0 ≤s≤Sh) 分別是單元軸向和周向的參數坐標.

圖1 薄壁桿件溫度單元及其單元局部坐標系Fig.1 thin-walled bar temperature element and its local coordinate system

忽略薄壁桿件內壁的輻射換熱,僅考慮其外壁向外輻射散熱.根據單位長度內“體內存儲的熱量 +體內傳導出去的熱量+表面輻射出去的熱量=表面吸收的熱量”可在單元局部坐標系中建立控制方程

式中,c是比熱容;ρ是密度;kξ和ks分別是桿件軸向和橫截面周向的導熱系數;αs和ε分別是桿件表面的熱流吸收率和輻射率;σ=5.67×10-8W/(m2·K4)是斯蒂芬-玻爾茲曼常數;Tbk是環境背景溫度.在軌環境中的太陽熱流非常均勻,不沿桿件軸線發生變化(除非其他結構的遮擋產生了陰影,但分析過程并不發生實質性的改變),所以有效熱流為

式中,n是外法線方向矢量;q=-qν 是入射熱流矢量(方向矢量 ν 從桿件表面向外為正方向);吸收熱流角 ψ(θ) 為 ν 和n之間的夾角,它沿薄壁桿件橫截面周向的變化規律取決于桿截面形狀.所以特定形狀薄壁桿件所吸收的有效熱流取決于各時刻入射熱流q(t)與該桿件橫截面各點的相對位置.

假設單元溫度沿軸向線性分布,則可采用兩節點Lagrange 形函數Ni(ξ) (i=1,2)在軸向插值.桿件橫截面內的溫度場可分解為平均溫度Tˉ 和攝動溫度兩部分,并將沿周向進行Fourier 級數展開

因為閉口截面的溫度場在周向呈周期性分布,而開口的斷面幾乎為絕熱邊界(即 ?T/?θ=0),所以上式中的周向插值函數Nn(θ) 為

根據式(4),每個單元節點包含1 個平均溫度和M個攝動溫度,其中攝動溫度反映了沿桿件橫截面的溫差.因為桿件橫截面的溫差(100量級)遠小于其平均溫度(102量級),所以方程(2)中

再注意到式(5)中的三角函數均滿足正交性條件[19-21],于是可以通過方程(2)和(3)得到將平均溫度和攝動溫度解耦的單元有限元方程.

單元方程經過組集后[19-21]得到的整體有限元方程仍然是解耦的,其中的平均溫度方程是非線性的

根據該方程解出平均溫度后,再求解如下線性方程

就可以得到攝動溫度.

基于上述理論框架,可以針對不同截面形狀的薄壁桿件開發出具體的單元[32-34],甚至可以考慮入射熱流被部分遮擋[33]等實際工程狀況.在某些特殊的熱流入射情況下,該方法還可以推廣至不等壁厚的薄壁桿[34].

1.2 熱特征時間

特定變形模式的熱誘發振動的劇烈程度和其對應溫度隨時間變化的快慢相關.通常沿桿件軸線方向的入射熱流分布是比較均勻的,所以只需要分析桿件橫截面內的平均溫度和攝動溫度隨時間變化的特性就能揭示熱誘發振動的機理.忽略方程(2)中沿桿件軸線方向的變化項后,再注意到式(6),可以得到解耦的平均溫度和攝動溫度的齊次方程

上述分析表明1.1 節中利用Fourier 級數正交性得到解耦的有限單元方程的方法不僅是一種數學技巧,而且符合熱傳導問題的物理特性.這種方法避免了同時求解慢變量(平均溫度) 和快變量(攝動溫度),既保證了求解精度又提高了求解效率.具體來說,方程(7) 僅包含了少量的慢變量成分,是“濃縮”后的純非線性部分.該方程可以用Galerkin 兩點差分格式進行時間積分,在每一時間步內用Newton-Raphson 法求解.而方程(8)雖然包含了大量的快變量成分,但由于已經從方程(7)中解出,所以該方程是線性的,可以高效求解.

1.3 溫度單元驗證算例[33]

算例1如圖2(a)所示,一根兩端絕熱、長度為100 mm 的閉口梯形薄壁桿受到被電池片部分遮擋的太陽熱流q(1)和地球熱流q(2)的照射并向外層空間輻射散熱.在初始溫度為290 K,環境溫度為0 K 的情況下,分別用ANSYS 軟件的四節點薄殼溫度單元和Fourier 桿件溫度單元計算其溫度場.在桿件橫截面周向均勻劃分了21 個薄殼單元,而Fourier 桿件單元中取M=6.如圖2(b)和圖2(c)所示,在此復雜情況下Fourier 桿件單元和薄殼單元的計算結果非常吻合.另外,同一截面上薄殼有限元模型一共有21 個溫度自由度,而Fourier 桿件有限元模型只需7 個廣義溫度自由度,從而減少了計算工作量.

圖2 薄壁梯形桿的溫度計算結果Fig.2 Temperature of a thin-walled trapezium cross-section bar

算例2如圖3(a)所示,將文獻[4]中哈勃太空望遠鏡的伸展桿由閉口圓管改為開口薄壁圓管,它受太陽熱流q的照射并向外層空間輻射散熱.初始溫度為290 K,環境溫度為0 K,分別用ANSYS 軟件的薄殼溫度單元(沿截面周向劃分24 個單元)和Fourier桿件溫度單元(M=3)計算其溫度場.圖3(b)和圖3(c)表明Fourier 有限元的結果有足夠的精度.圖3(d)給出各諧攝動溫度的時變歷程,和圖3(b)相比,可以明顯地看出攝動溫度的變化速度遠大于平均溫度的變化速度,驗證了本文1.2 節的分析結果.

圖3(c)顯示當入射熱流方向與開口薄壁桿件的幾何對稱面不重合時,沿桿件橫截面周向溫度分布不對稱,因此不僅造成彎曲變形,還引起非對稱的截面翹曲變形,從而產生熱雙力矩使桿件發生扭轉.

圖3 薄壁C 形桿的溫度計算結果Fig.3 Temperature of a thin-walled C-shaped cross-section bar

2 熱-動力學耦合有限元分析

Fourier 溫度桿件單元可以和常規的結構動力學桿件單元協調工作,從而便于分析熱誘發振動問題.熱-動力學耦合效應分為弱耦合和強耦合兩種情況.弱耦合只關心瞬態溫度場所產生的結構動態熱變形,而強耦合還需考慮結構熱變形改變了它所吸收的有效熱流后對溫度場造成的影響.分析熱顫振問題時,必須考慮熱-動力學強耦合效應[10].

柔性可展開附件的結構動力學分析應該考慮幾何大變形的影響.對于剛性太陽翼等比較剛硬的結構,則完全可以采用線性方法進行分析.

2.1 熱-動力學耦合有限元分析的一般框架

最一般的情況,需同時考慮輻射換熱、幾何非線性和熱-動力學強耦合三種非線性問題.此時應采用迭代法逐步求解[20-21,33].如圖4 示,設t時刻的所有場變量都為已知,則可在全局坐標系(X,Y,Z)中求解t+Δt時刻的節點溫度t+ΔtTg和節點位移t+Δtug.

圖4 桿單元的變形過程Fig.4 Deformation of the bar element

首先用第1 節的方法求解t+ΔtTg.為表述方便,可認為是求解方程(7)和(8)的合并形式

上式中的熱流載荷向量是先在t時刻單元局部坐標系(tx,ty,tz)中根據式(3)計算,然后再轉換到全局坐標系(X,Y,Z)中的.

在假設薄壁桿件的橫截面是剛周邊的前提下,定義在單元局部坐標系中的桿件外表面法向矢量n沿桿件橫截面周向的分布只和桿件橫截面形狀有關[37-38],不隨桿件變形而改變.但是,由于太陽熱流方向是恒定的,定義在單元局部坐標系中的入射熱流方向矢量tν 將隨著桿件的變形tug而改變(具體可見2.2 節).

根據t+ΔtTg可計算結構所受的熱載荷向量t+ΔtF(包括了由平均溫度產生的熱軸力,攝動溫度產生的熱彎矩和熱雙力矩),然后通過如下更新拉格朗日格式的非線性動力學方程求解t+Δtug

其中,t f是內力向量;tM和t D分別是質量矩陣和阻尼矩陣.由于熱誘發振動是典型的小應變大轉動問題,所以就是常規的線彈性剛度矩陣.而幾何剛度矩陣包含了初應力的貢獻,也考慮了溫度變化的影響.方程(11)可用Newmark 法進行時間積分,并在每一時間步內用Newton-Raphson 法迭代求解.

結構變形后其所受熱流會發生變化,從而導致溫度改變,并進一步影響熱變形.這種熱-變形的強耦合效應在熱傳導方程(10)中體現為t+ΔtQ和位移tug相關,動力學方程(11)中t+ΔtF和溫度t+ΔtTg相關.為保證解的正確性,需要迭代求解方程(10)和(11),直至t+ΔtQ的變化量小于預設值.迭代過程中,這些方程中的矩陣和向量都需要根據構型的變化進行更新.

在求解過程中還需要注意以下幾個問題.

(1) Newton-Raphson 法的切線預測步驟所用的斜率是比較靈活的,所以在KNL中只考慮桿件軸力引起的非線性效應即已體現關鍵物理特征[33].

(2) Newton-Raphson 法的路徑修正步驟需要正確地計算不平衡力.若直接用積分方法計算內力向量,會因構型轉換而非常繁瑣.而根據剛體轉動準則[39],考慮到Cauchy 應力在隨體坐標系下不會因為剛體轉動而改變,就可以直接將k-1 構型單元局部坐標系中的內力k-1fe(而不需要轉換到k構型)加上表達在k構型局部坐標系中的增量內力 Δfe,得到k構型單元局部坐標系中的內力

(3) 對于開口薄壁桿需要考慮截面翹曲位移所產生的約束扭轉效應,以及桿件橫截面形心和扭心之間的坐標轉換關系[33-34,40-41].

(4) 應采用2.2 節介紹的Rodrigues 公式來描述桿單元的空間大轉動關系,以保證求解精度.

2.2 大轉動情況下的構型轉換關系[33]

方程(10)和(11)中的矩陣和向量是首先在某個構型的單元局部坐標系中建立后,再轉換到全局坐標系(X,Y,Z)中的.這些坐標轉換會直接影響求解精度,特別是熱-變形強耦合分析的準確性.

如圖4 所示,以位移為例,因為初始構型單元局部坐標系(0x,0y,0z)是已知的,所以可以用常規的坐標轉換矩陣將初始構型中的向量0ue表達為全局坐標系(X,Y,Z)中的向量0ug

將初始構型的單元局部坐標系(0x,0y,0z)進行剛體轉動后就可以得到t時刻構型的單元局部坐標系(tx,ty,tz),從而可以用相應的正交矩陣建立0ue和tue之間的關系.但是對于熱誘發振動這種小應變大轉動問題,若用Euler 轉動公式來計算(tx,ty,tz)不但失去了保內積的性質,而且所得結果還依賴于繞坐標軸的轉動順序,因此必須使用精確的正交矩陣進行計算.如圖5 所示,此時認為(0x,0y,0z)是先繞正交矩陣的軸方向矢量teβ旋轉tβ 角度到達(tx,坐標系后,再繞tx軸旋轉tγ 角度到達t時刻單元局部坐標系(tx,ty,tz).相應的轉換關系為

圖5 桿單元的構型轉換Fig.5 Configuration change of the bar element

式(15)是精確的,和繞坐標軸的轉動順序無關.當tβ 很小時,式(15)就退化為Euler 轉動公式.根據式(13)和式(14),就可以實現在全局坐標系和任意時刻單元局部坐標系之間的轉換.

2.3 熱-動力學強耦合算例[21,34]

采用上述全耦合方法分析圖6 所示哈勃太空望遠鏡的熱顫振事故.文獻[4]所建立的傳熱學和力學模型是閉口薄壁圓桿,只能模擬熱誘發彎曲振動現象,無法回答“熱為什么會引起扭轉振動”的問題.實際上該太陽翼的Bi-Stem 支承梁是開口薄壁桿件,在軌展開后其開口方向是隨機的.為了更加接近這種實際情況,在有限元模型[24,40-41]中采用開口C 型圓桿來模擬支承梁,并假設左右梁的開口方向有θ=20°的偏差.太陽翼中部的柔性太陽毯采用等效的閉口薄壁圓桿單元網格模擬.模型中還考慮了在軌展開后太陽毯受到的29.5 N 的預張力,以及左右梁受到的14.75 N 預壓力.模擬時取太陽熱流為1350 W/m2,沿 ψ=80°角入射.

圖6 哈勃太空望遠鏡太陽翼模型Fig.6 The model of the HST solar array

從圖7 可以看出,熱誘發振動包含了彎曲(0.097 Hz)與扭轉(0.004 0 Hz)兩種不同模式的振動.大約在1200 s 時太陽翼開始失穩,此時梁的壓縮軸力達到了壓桿失穩的臨界值.圖7(c)顯示的失穩后的彎扭耦合變形圖和實際照片非常相似.由此證明了本文方法的有效性.圖7(a)還說明,雖然采用線性分析可以模擬出彎扭耦合熱誘發振動,但是無法得到熱顫振和熱屈曲的正確結果.可見對于柔性大型空間結構,必須采用非線性分析.

圖7 哈勃太空望遠鏡太陽翼的模擬結果Fig.7 The simulation results of the HST solar array

3 產生熱誘發振動的必要條件

式(1)只適用于評價突加熱流作用下純彎曲梁的熱誘發振動劇烈程度.而實際工程結構的振動模式是非常復雜的,所受熱流也不是理想的階躍函數.故需討論一般情況下產生熱誘發振動的必要條件.

3.1 理論分析[37]

為了能得到類似于式(1)的解析定量結果,本節和Boley[7-8]原始論文一樣,只采用線性結構模型且不考慮熱-結構強耦合效應.為得到一般的結果,采用了有限元模型進行分析,此時方程(13)退化為

式中的剛度矩陣K可以包含預應力.熱載荷FT包含了熱軸力、熱彎矩和熱雙力矩.其中熱軸力和桿件橫截面平均溫度有關,是慢變量;其他和桿件橫截面的攝動溫度有關,是快變量,可能導致熱誘發振動.為評估產生熱誘發振動的必要條件及劇烈程度,用模態疊加法求解由攝動溫度引起的位移[37].

記為產生的熱誘發振動位移,它在振動模態空間表達為分別是系統的第i階振動模態和模態坐標.將熱載荷轉換至振動模態空間,于是方程(16)可以解耦為I個獨立的方程

式中,ωi和 ?i分別為系統的第i階(i=1,2,···,I)振動圓頻率和阻尼比,而 Θix(t) 可以表示為

于是可知在突加恒定熱流并忽略阻尼的條件下,第i階模態熱誘發振動的劇烈程度為

上式說明第i階模態熱誘發振動的劇烈程度不但和 ωi有關,而且和各階熱特征時間 τj都有關.若只考慮第一階熱特征模態(J=1),則式(21)退化為式(1).式(21)對具有復雜模態的工程結構仍成立,而Boley[6-7]所給出的熱誘發梁彎曲振動只是一特例.

式(21)成立的條件也是非常理想的.實際結構中會存在阻尼,可能是非線性的,所受熱流也不是理想的階躍函數.此時式(21)的估計往往是保守的.

3.2 算例說明

算例1文獻[4]將哈勃太空望遠鏡太陽翼的Bi-Stem 梁(見圖6)簡化為閉口薄壁圓管,給出了其熱誘發振動的理論估計.文獻[43]改變帆板伸展梁的厚度和傳熱參數,用本文數值解得到不同的Boley參數與帆板端部最大位移占準靜態位移之百分比的關系(見圖8),和理論預測非常吻合.

圖8 桿端部的熱誘發振動劇烈程度隨B 的變化Fig.8 The bar tip TIV intensity changes with B

算例2圖9 是一個典型的半剛性太陽翼力學模型,主要由受18.9 N 預張力的張緊網、受預壓力的支承梁、頂板以及邊框等組成.整個結構由10 塊太陽能帆板展開而成,每塊帆板分為4 個700 mm ×750 mm 框.帆板框由薄壁正方形管構成,支承梁為薄壁圓管.帆板框和支承梁固結于衛星艙體上.有限元分析時將電池片質量計入由凱夫拉纖維構成的張緊網纖維上,略去艙體運動對帆板的影響并假設二者之間絕熱.從0 時刻開始,結構受到垂直于帆板平面的太陽熱流照射,并向太空輻射散熱.太陽翼各構件的尺寸和材料見表1.

圖9 突加熱流作用下的半剛性太陽翼Fig.9 A semi-stiff solar panel subjected to suddenly applied heat flux

表1 太陽能帆板各構件的尺寸及其材料Table 1 Dimension and material of the solar panel

文獻[43]采用桿件Fourier 有限元對該太陽翼模型進行離散,總計682 個節點,1052 個單元.其中支承梁采用薄壁圓管單元,張緊網繃弦采用實心圓桿單元(其材料參數已經等效了電池片的質量和熱容量),帆板框和頂板為薄壁矩形管單元.分別用直接積分法和基于Lanczos 方法的模態疊加法計算溫度場,所得攝動溫度吻合得很好(見圖10).用Lanczos 方法得到結構1 邊框的熱特征時間為10.47 s,而對攝動溫度時程曲線以指數函數擬合得到的熱特征時間約為10.5 s,二者一致.

圖10 帆板框(結構1)和支承桿的攝動溫度時間歷程Fig.10 Change of perturbation temperatures of the frame and boom

圖11 是采用弱耦合方法得到的A點(見圖9)熱誘發振動響應.其中的結構1 和結構2 可見表1.表2 顯示,對于這個比較復雜的結構實際的熱誘發振動劇烈程度比式(1)的估計低.這是因為該結構的第一階模態不是理想的純彎曲梁形式(帆板是板振動模式,和支承梁的彎曲模式不同),熱特征時間也不能完全由邊框代表.

圖11 突加熱流作用下太陽翼端部的法向位移Fig.11 Tip displacement of the solar panel subjected to suddenly applied heat flux

表2 A 點的主要特征量 (t=2000 s)Table 2 Main features of point A at t=2000 s

算例3實際的入射熱流總是經過一段時間后才逐漸穩定的,不可能是階躍函數.因此經典的熱特征時間總是小于實際溫差趨于穩定的時間,熱誘發振動的劇烈程度也會小于式(1)或式(21)的估計值.比如,在軌航天器在進出地球陰影時,總會經過一段半影區,軌道越高半影區越寬,熱流變化越慢,越不容易產生熱誘發振動[3].另外,柔性結構由于發生了大轉動,吸收的熱流會比只發生小轉動的剛性結構少一些,熱誘發振動也會小一些(沒發生熱顫振時).文獻[28]根據1.5 m 長和2 m 長之柔性懸臂梁的熱誘發振動試驗的測量數據,進一步討論了這些因素,并根據溫差變化情況擬合出等效的熱特征時間τeff和相應的等效Boley 參數Beff=ω1τeff.從圖12 的比較結果來看,經典理論預測的熱誘發振動最劇烈,等效理論的預測結果和實測結果接近.

圖12 經典理論、等效和實測振動比的比較[28]Fig.12 Comparison among the original,effective and actual ratio of dmax/ dst[28]

4 熱誘發振動的穩定性

熱誘發振動發生不穩定的物理機制是結構振動變形使結構所吸收的有效熱流qˉ 因結構振動變形而震蕩變化(見式(3)),從而產生熱-結構系統的運動穩定性問題.穩定性是系統略微偏離平衡狀態后表現出的性質,所以應該在結構變形之后的構型上進行討論.另外,結構發生失穩說明產生了不唯一的解,意味著這一定是非線性問題.但根據里雅普諾夫關于運動穩定性問題的一次近似判別準則[44],可用系統方程在變形之后構型的一階Taylor 展開項的性質來判別是否失穩.本節先借助懸臂桿件熱顫振這個簡單問題說明研究熱顫振的基本方法,然后再將它推廣至復雜的工程結構.

4.1 懸臂梁熱顫振準則

以圖13 所示模型來解釋引言中介紹的熱顫振準則.以半徑為R的閉口薄壁圓管為例,根據第1 節的理論可得到解耦的平均溫度和一階攝動溫度的控制方程

圖13 懸臂梁受突加熱流作用Fig.13 A cantilever beam subjected to suddenly applied thermal flux

對無阻尼的懸臂桿,令z方向位移w的插值形式為w(x,t)=N(x)W(t),方程(16)變為

采用同樣的思路,也可以得到太陽帆的熱顫振條件[45],甚至在考慮開口方向的影響后還可以進一步得到開口桿件的熱顫振準則[34,46].

4.2 熱顫振的能量解釋

發生熱誘發振動時,入射熱流的方向ν并未改變,所以這是一種典型的自激振動.當系統從外界吸收的熱量超過發散的熱量時,就會發生熱顫振.因此,也可以從能量的角度推導出懸臂梁的熱顫振準則[29,47].

圍繞著穩態平衡位置,梁的轉角是振蕩的θ(b)(t)≈相應的法向有效熱流為

將式(26)代入式(23)并注意到ks/R2?4εσTˉ3/ts,可得

一個周期T=2π/ω 內,熱彎矩做的功為

4.3 復雜結構的熱顫振分析[38]

對復雜工程結構可先計算其低階固有頻率和對應的熱特征時間,按照第3 節“發生熱誘發振動的必要條件”,初步判斷是否會發生熱誘發振動.

對于柔性結構,需按第2 節,計算結構振動隨時間的變化,判斷會否發生熱誘發振動與熱顫振.

若結構剛度較大,可按有限元方程直接分析熱誘發振動的穩定性區間[34].將式(17)寫成矩陣形式

并在式(8)中考慮振動耦合項BIU的影響

式中,J為相對于穩態的一階近似;ν為入射熱流的方向矢量; ? 為阻尼比向量.

根據常微分方程理論,當矩陣L(ν,?)的最大特征值的實部α＜ 0 時,系統是穩定的;當α＞ 0,則系統不穩定;當α=0,則系統處于臨界狀態.

現以一個典型的環形天線(見圖14)在突加太陽熱流作用下的穩定性分析[38,48]來說明本文方法的實際應用.環形天線由主桿、圈桿、豎桿和上下反射網構成,其截面尺寸以及材料見表3.天線各構件都采用薄壁圓管單元模擬,總計748 個單元.假設天線通過主梁固結于衛星艙體,并忽略與衛星之間的傳熱,只研究環形天線本身的熱-動力學問題.

表3 環狀天線的單元類型,尺寸與材料Table 3 Dimension and material of the circular antenna

該天線系統的穩定性取決于三個參數: 入射熱流矢量的方位角φ1和φ2和系統阻尼系數?.根據式(35)矩陣的最大特征值實部α(φ1,φ2,?)=0 的條件,可解得臨界曲面,它將三維參數空間劃分為穩定區和不穩定區(見圖15).

利用圖15 所給分區,分別在穩定區、臨界面和不穩定區計算環形天線的熱誘發振動響應,可進一步驗證本節運動穩定性分析的可靠性.圖16 是在給定熱流方向時阻尼比和α的關系.可見阻尼小于臨界值?cr后系統就不穩定,此結論也可通過圖17 中A點(見圖14)的Z向振動情況得到印證.

圖14 突加太陽熱流作用下的環狀天線Fig.14 Circular antenna under the incident heat flux

圖15 環狀天線的穩定性邊界Fig.15 Stability boundary of the circular antenna

圖16 α 的根軌跡圖Fig.16 The locus of α

圖17 A 點Z 向振動Fig.17 The vibration in Z direction at point A

5 其他實際問題

前文已經將熱誘發振動和熱顫振的機理、特性和分析方法做了詳細的介紹,其關鍵在于理解慢變量和快變量性質,并采取相應的應對措施.本節將進一步討論分析復雜工程結構時需要注意的問題.

工程結構往往由很多不同種類的構件組成.對于不同截面形式的桿狀構件,需要用不同的Fourier桿件溫度單元.另外,還有一些構件不是桿狀的,比如太陽翼的電池板.為了能和Fourier 桿單元配合使用,可以采用等效的梁格模型[34,49],即用Fourier 桿單元網格來等效電池板.此時不但要求所有力學特性(比如頻率)等效,還要求熱學特性(比如熱特征時間)等效.

現代航天器的可展開結構不斷向大型化和柔性化方向發展.當它們的轉動慣量和航天器艙體相比不再是一個小量時,就必須考慮和艙體的耦合動力學問題[38,49].此時可以采用強耦合方式進行分析.但穩定性分析的結果表明,空間結構和航天器艙體的轉動慣量比越大,系統越容易失穩,此時需要再考慮和姿態控制系統的協調關系[44].

工程中還特別關心怎樣避免產生熱誘發振動,或者怎樣抑制熱誘發振動.由于熱誘發振動的頻率很低(通常在0.1 Hz 左右甚至更低),被動隔振方法效果不明顯,主動控制不但難度很大而且還要付出高能耗和可靠性方面的代價[50-51].所以采用優化設計方法,從源頭上盡量避免發生熱誘發振動,更便于工程應用.最簡單的措施就是通過增加反射涂層來減小輸入熱流.另外,也可以根據特定的需要做更復雜的優化設計[52-53].當然,最根本的解決辦法是采用零膨脹材料來制造空間結構[24].但是在軌航天器的表面溫度的變化范圍在正負一百多度之間,怎樣保證在這么大的溫度范圍內材料都有穩定的零膨脹性能,也是極具工程挑戰性的任務.

由于熱誘發振動產生的條件非?？量?所以驗證試驗的難度非常大.理想情況應該進行在軌試驗[24],但要付出巨大的代價.所以通常還是開展地面試驗,即便如此其難度和花費都是很高的[26-29].這種試驗需要在熱真空罐中開展,用紅外燈陣模擬太陽熱流,并最好用抽取遮陽板的方式來模擬從陰影到光照區的加熱過程.此外在試驗過程中要特別注意以下幾點.

(1) 和空間太陽熱流相比,在有限距離內紅外燈陣產生的熱流方向不是平行的,其分布密度還隨燈與試件的距離而變化.因此需要考慮這些因素對試驗結果的影響[3,27-29].

(2) 直射在熱電偶上的熱流會影響對試件表面溫度測量的準確性,從而降低了薄壁構件截面溫差測量結果的精度.因此需要精心設計熱電偶的安裝和溫度補償措施.

(3) 瞬態熱流很難直接測量.通常是用測量溫度根據熱傳導方程反算出熱流隨時間變化的曲線.但是由于熱電偶測量的延遲效應,反演的熱流變化速度會比真實情況慢.此時可通過熱誘發振動的幅值對輸入熱流進行修正[28].

(4) 重力對柔性結構地面試驗結果的影響很大.但由于要求試驗必須在真空中進行,因此不能采用氣浮裝置消除重力.通?？刹捎脗认驊覓旎蜇Q直懸掛的方式來部分解決此問題.但是豎直懸掛試件的振動幅值會因為重力產生的恢復力而遠小于真空中的幅值.側向懸掛受重力的影響會小一些,但是因為懸掛裝置不可能無限長,所以也會產生額外的約束,導致試驗結果失真[28].

總之,由于熱環境和重力等各種條件的影響,地面試驗結果不能完全和在軌的真實結果吻合.在無法開展在軌試驗的情況下,現在只能采用數值模擬和地面試驗相結合的方法來預測在軌的熱誘發振動響應[28].即,用地面試驗的結果修正有限元模型,然后用有限元方法分析在軌真實環境下(比如失重狀態)的結構響應.因此有限元模型和試驗測試的精度就至關重要.

材料參數的準確性對數值模擬結果有很大的影響[54].但是一些熱學參數測量的分散性比較大(特別是復合材料),而且力熱參數通常是隨溫度變化的.采用隨機有限元方法可以考慮材料參數的不確定性,但計算工作量很大,此時可以采用響應面方法來提高分析效率[53].若要在模擬過程中考慮參數隨溫度的變化,則需要增加額外的計算工作量.比如在隱式算法中增加額外的迭代步驟,或者直接采用時間步長較短(保證數值穩定性)的顯式算法.不過,考慮到材料參數隨溫度變化的因素也是一個慢變量,而且熱誘發振動往往是在進出地球陰影的較短時間內發生的(此時平均溫度還沒有發生太大的變化),那么也可以取所關心的溫度范圍內的參數平均值進行分析.但是材料參數的不確定性所造成的誤差是否小于空間望遠鏡等高精度設備的要求,還是需要進一步討論的.

6 結論

經過幾十年的努力,人們雖然已經理解了熱誘發振動現象的機理,但是在工程應用時還是面臨著諸多挑戰[6,23].特別是在研制空間望遠鏡等對指向性要求特別苛刻的航天器時,應該綜合考慮光-熱-結構-控制系統的一體化耦合分析[54].此時,將多個現有軟件系統集成的方式已經不能滿足這種高精度的要求,而應該發展專門的一體化數值模擬方法和與之配套的地面試驗和驗證技術.這正是未來需要努力的方向.

致謝

20 多年來,清華大學先后有多位老師與博士、碩士研究生參加了本課題組的工作或給予了寶貴的建議,他們是胡寧、任革學、岑章志、曹陽、丁勇、程樂錦、姚海民、黃彥文、李軍、張曉東、張逸凡、唐羽燁、段進、李偉、范立佳、張軍徽、樊孝清、袁小德和金鄧格.他們的創造性研究和辛勤勞動,為課題的順利完成做出了重要貢獻,作者謹在此致以衷心的感謝.