熱-動力學耦合多體系統建模與降階1)

田青龍 於祖慶 蘭 朋,,,2) 陸念力

* (哈爾濱工業大學機電工程學院,哈爾濱 150001)

? (河海大學機電工程學院,江蘇常州 213002)

** (西安建筑科技大學機電工程學院,西安 710055)

?? (西安建筑科技大學西部綠色建筑國家重點實驗室,西安 710055)

引言

太陽能電池板和拋物面天線等典型的航天器板結構具有重量輕、尺寸大的特點.此類航天器的動力學分析與控制問題越來越多地需要考慮到空間環境溫度的影響.由太陽輻射導致的極端熱載荷激勵引起位移場與溫度場的強耦合響應,進而導致熱致振動或熱跳變[1-4].因此,剛柔熱耦合系統的精確建模對航天器的安全至關重要.

早期的熱動力學耦合研究較多使用浮動坐標方法(FFRF)[5-6],因而其僅限于柔性體小變形問題.絕對節點坐標方法(ANCF)則采用了全局坐標系下定義的位置與梯度向量作為節點坐標,被廣泛用于大范圍運動和大變形的多柔體系統動力學的研究[7].Wu等[8]在假設梁截面的溫度分布為線性的情況下,研究了瞬態熱載荷對細長梁結構的影響.Abbas 等[9]使用變厚度ANCF 板單元研究了熱載荷下的氣動熱彈性動力學.?epon 等[10]建立了基于二維ANCF 梁單元的雙金屬片模型.在上述研究中,熱載荷均以解析解的形式給出.Li 等[11]建立了熱柔耦合的大柔性復合太陽能電池板平面航天器模型,并對其進行了動力學和熱耦合分析.于香杰等[12]結合負應變率控制法建立了變截面梁構件在空間熱載荷作用下的壓電抑振抑變結構.邢曉峰等[13]研究了剛-柔耦合、耦合熱彈性、耦合熱-結構三重耦合效應對航天器結構熱致振動響應及穩定性的影響.Shen 等[14-15]研究了薄壁管和多層材料太陽能電池板的熱致振動.溫度場采用線性插值網格進行離散且傳熱方程與動力學方程同步求解.這一階段的研究均采用兩種不同的離散格式分別對柔性體位移場和溫度場進行插值,會導致數值求解中的網格錯配和映射誤差問題[16].為此,Cui 等[17-18]將ANCF 單元形狀函數應用于位移場和溫度場,自然地將傳熱和連續介質力學統一為一個單元網格.這一思想也將在本研究中得到應用.

隨著ANCF 計算規模增大,數值效率低的問題日漸突出,模型降階成為了重要的研究方向.在文獻[19-20]的研究中,本征正交分解(POD)方法被用于ANCF 單元的模型降階.然而,降階模型必須通過實驗數據庫或計算全階模型仿真結果來獲得.模態綜合法(CMS) 可以避免這一缺點.Kobayashi 等[21]在軸向變形較小的前提下,使用固定界面模態綜合法(也稱Craig-Bampton 方法)來縮減ANCF 平面梁單元建模的多體系統.Sun 等[22-23]使用自由界面模態綜合法來提高帶滑動鉸的柔性梁的計算效率.該方法同樣采用了相同的軸向小變形的假設.Tang 等[24]將整個柔性體運動過程劃分為若干子區間并在每個子區間內對系統方程進行線性化處理,從而使得切線剛度陣定常,基于模態截斷的降階方法得以應用.此外作者也采用了Craig-Bampton 方法對ANCF 梁單元描述的多柔體系統進行子結構劃分,稱之為CB-ANCF 方法.在此基礎上,作者將此方法推廣到了具有阻尼的多柔體系統,以減少黏彈性多柔體系統的自由度[25].基于文獻[24-25] 的研究,Tian等[26]使用自由界面模態綜合法來減少ANCF 梁單元描述的多柔體系統的自由度,即F-ANCF 方法.

與動力學方程不同,傳熱方程為一階偏微分方程.很多研究表明,傳統的動力學降階方法如模態疊加法[27]、模型識別法[28]、POD 法[29]等依然適用于傳熱方程的降階.然而,對于熱柔耦合系統的整體降階研究卻較少見于文獻.Nachtergaele 等[30]提出了一種弱耦合投影基用于熱柔耦合系統的降階,這是對Craig-Bampton 方法的擴展,溫度場隱含在投影基中.Yamashita 等[31]也使用該方法對基于FFRF 建模的熱柔耦合系統進行了降階研究.

本文提出了一種剛柔熱耦合系統模型降階的方法.引入溫度梯度,采用高階插值對溫度場進行離散,使位移場和溫度場采用統一的ANCF 單元形函數進行插值離散,避免數值求解中網格錯配和誤差映射的問題.ANCF 參考節點用于描述中心剛體.采用泰勒展開對動力學方程和傳熱方程進行線性化處理,得到常值切線剛度陣.引入改進的模態綜合法分別對動力學方程和傳熱方程進行自由度縮減.最后,以剛柔熱耦合大口徑拋物面天線為例驗證了該方法的有效性.

1 熱-動力學耦合建模

本文擬分別采用熱耦合的ANCF 薄板[18]與三維梁單元[17]分別構造星載拋物面天線及其肋梁模型.而衛星本體視為剛體,由ANCF 參考節點[32]描述.由此可將衛星剛柔熱耦合系統納入ANCF 框架內描述和建模.由于采用了全局坐標系下定義的位置與梯度向量,單元內部任意點的空間位置可以定義為r=Se,S為單元形函數矩陣,e為單元節點向量.ANCF 薄板單元基于Kirchhoff-Love 板理論,其節點向量包含中面四個角點的位置和平面內梯度向量,即其中上標P 代指板單元,i=1,2,3,4 為節點號.ANCF 三維梁單元節點為中軸線兩端點,每個節點使用全局位置向量和沿三個方向的梯度向量作為節點坐標,可以表示為:eB=上標B 代指梁單元,j=1,2.

ANCF 參考節點是一種具有完備梯度向量但不屬于任何單元的孤立節點.它通過連續性條件與柔性體網格相連,并通過非線性約束的形式保證其梯度向量模長始終為1 并且正交.因此它可以在ANCF的表達范式下高效地描述剛體的動力學行為.系統動力學方程的形式不因其的加入而發生改變.現有的研究表明,無論柔性體離散使用全參數單元還是梯度不完備單元,均不影響ANCF 參考單元描述剛體運動的正確性.因此,考慮熱應力的系統動力學方程為

其中,M為系統質量矩陣,Qe是彈性力.薄板單元和梁單元的彈性力分別在文獻[33-34]中給出.QH是對應于熱應變的廣義力,表達式在附錄中給出.Ce為約束方程,?Ce/?e則為約束方程對廣義坐標e的偏導數,λ 是拉格朗日乘子,Qext代表廣義外力.

由于連續體內溫度場采用了與運動學描述相同的離散格式,單元節點處的溫度及其梯度被用作傳熱方程的廣義坐標,如圖1 所示.薄板單元的溫度廣義坐標為其中Ti(i=1,2,3,4)為節點i的溫度,Txi和Tyi分別表示節點i在x和y方向上的溫度梯度.同理,三維梁單元的溫度廣義坐標可以寫成

圖1 ANCF 熱耦合單元Fig.1 ANCF thermal coupling element

與位置插值同理,單元內任意點的溫度可以寫為

其中,ST是溫度場的單元形狀函數,其元素與描述位移場形狀函數相同,eT表示單元的節點廣義溫度坐標向量.同樣地,上角標P 和B 分別指代薄板單元和梁單元.

在太空中,影響航天器溫度變化的主要原因是太陽熱流輸入和表面自熱輻射,因此在考慮熱的邊界條件的情況下,傳熱方程的一般形式如下

其中,DT是熱容矩陣,Kc是熱傳導矩陣,Kr是表面輻射矩陣[17-18].eT∞代表熱輻射環境溫度,在太空中自熱輻射環境溫度近似為0 K.QT是外熱載荷,與熱邊界條件有關.當受到太陽熱流照射時,外熱載荷通常與柔性體當前構型和太陽熱流入射角度有關:QT=∫H·ndS,H=H0nT,n為受熱面法線方向,nT為熱流入射方向,H0為太陽熱流密度.CT傳熱方程的約束方程. ?CT/?eT是傳熱約束方程對廣義溫度坐標eT的偏導數,λT是傳熱約束的拉格朗日乘子.

2 動力學和傳熱方程的線性化

基于ANCF 建立的柔性體動力學方程具有常值質量陣的優勢,但其彈性力卻是節點坐標向量的高度非線性函數,切線剛度陣亦然.這就導致基于模態截斷的降階方法無法直接應用.為此,將系統的整個運動過程劃分為多個區間,認為柔性體在每個區間內部的動力學行為是準靜態的.利用一階泰勒展開將基于初始構型定義的全量平衡方程變換為基于當前區間起始構型的增量平衡形式.每個區間內的切線剛度矩陣可視為定常,進而采用改進的模態綜合法對系統進行降階.

根據泰勒展開式,每個準靜態區間的初始構型處的動力學方程可展開為

其中,e0是每個準靜態區間初始構型處的廣義坐標.同理,eT0是初始廣義溫度坐標.JQe是彈性力雅可比矩陣(表達式見文獻[33-34]).JQH為熱應變廣義力的雅可比矩陣,其表達式在附錄中給出.約束方程也可以線性化為

將式(4)和式(5)代入式(1)中,線性化的動力學方程形式可以寫成

其中,K0是廣義剛度矩陣,F0是線性化的廣義外力,表達式為

用同樣的方法,傳熱方程可以展開如下

傳熱約束方程可展開為

將式(7)和式(8)代入式(3)中,線性化的傳熱方程形式可以寫成

式中,FT0是線性化的廣義熱載荷.

3 改進的模態綜合法

根據子結構界面約束的不同,模態綜合方法可分為固定界面模態綜合法[24]和自由界面模態綜合法[35].如圖2 所示,多體系統ABCD被連接節點EF分為兩個子結構.對于固定界面模態綜合法,子結構之間連接的自由度被全部保留.然后固定每個子結構的界面自由度,并分別對子結構進行模態分析.最后,子結構之間通過共享邊界節點的形式連接起來,如圖2(b)所示.而自由界面模態綜合法則無需保留連接界面自由度,即每個子結構的模態都包含了剛體位移信息.結構的組裝通過界面位移和界面力相等的雙協調條件實現,如圖2(c)所示.而本文在自由界面模態綜合法的基礎上提出了一種改進的模態綜合法.通過施加界面位置和梯度的約束方程保證子結構間的連接精度,如圖2(d)所示.

圖2 系統子結構劃分Fig.2 The mesh of the system and substructures

在動力學方程部分,子結構j的動力學方程齊次形式如下

其中,[]j表示子結構j.通過模態分析獲得當前子結構的特征值和特征向量

其中,ωn是第n個固有頻率,[Φ]j為子結構振型,是子結構j的自由度數.

子結構j的廣義坐標增量可以表示為

以兩個子結構組成的系統為例,其廣義坐標可以表示為

顯然,每個子結構的坐標是獨立的,除了連接處的自由度.在本研究中,直接施加界面位置和梯度的約束方程,以確保子結構連接處的精度,即

將式(14)和式(15)代入式(6),系統動力學方程可以表示為

其中N=diag([Φk]1,[Φk]2)

傳熱方程是一階微分方程.子結構j的特征值和特征向量可以通過熱特征值分析獲得

式中,ωm為第m熱模態頻率,物理意義為熱衰減系數[27],[Ψ]j是對應于熱特征值的特征向量是子結構j的溫度坐標數.傳熱方程的特征值的求解在文獻[36]中給出.

子結構j的廣義溫度坐標增量可以表示為

直接施加界面連接處溫度和溫度梯度的約束方程

然后,將式(21)和式(22)代入式(9),系統傳熱方程可以表示為

4 求解策略

為準確模擬熱-動力學耦合效應對系統的影響,需要同步求解傳熱方程和動力學方程.本文采用廣義-α 積分來求解傳熱方程和動力學方程.圖3 給出了算法的流程圖.算法的數值參數 αf,αm,β ,γ 等的取值見參考文獻[37-38].

圖3 算法求解流程Fig.3 The overall flow-chat of the computation algorithm

在第2 節中,使用泰勒展開法將動力學方程和傳熱方程線性化,將非線性方程轉化為線性方程.需要指出的是,線性方程等式成立的前提是廣義位置坐標的增量應該很小.當第j個子結構的最大位移增量過大時,第j子結構應重新線性化.同樣,溫度的升高也會引起熱應變的變化,同樣會影響切線剛度陣,因此,當第j子結構的溫度增量過大時,第j子結構也應重新線性化.保證 δk和 δT足夠小,將累計誤差控制在一定范圍內.本文給出閾值的經驗計算公式

其中L是梁結構的長度或者薄板結構長度.通常,對于動力學方程,保留的自由度nred建議如下nred=ntal/3～4ntal/5.對于傳熱方程,保留的自由度建議如下

需要指出的是,t=0 時的初始模態速度和加速度均為0,任意子結構j第一次縮減保留的模態為前k階低階固有頻率對應的主模態.當 ‖Δr‖2＞δk或者 ‖ΔT‖2＞δT時,動力學方程和傳熱方程需要重新線性化,更新系統的切線剛度陣.之后,用新的切線剛度陣和質量陣去重新進行子結構模態分析,然后根據模態速度的絕對值大小排序,保留模態速度最大的前k階對應的主模態.更多細節可參見文獻[26].

5 數值算例

本節給出了四個數值算例來驗證所提出方法的準確性和有效性.為方便起見,完整自由度模型記為FOM-TANCF,降階模型記為ROM-TANCF.

5.1 半圓薄板的純導熱

在半圓形薄板上測試純導熱問題,一側受到熱流輸入,并考慮了表面自熱輻射.如圖4(a)所示,將半圓板劃分為四個ANCF 熱薄板單元.圖4(b)是參考文獻[14]中給出的一維線性傳熱單元,將半圓劃分十個單元.總仿真時間設置為80 s,初始溫度為290 K,熱流為1350 W/m2.其他參數見表1.

圖4 半圓板傳熱模型Fig.4 The heat conduction of a cylinder thin plate

表1 半圓板模型參數Table 1 Parameters of the semicircular thin plate

由于半圓形結構簡單,因此不進行子結構劃分.圖5 給出了FOM-TANCF、線性傳統有限元(FEM)和解析解的結果.可以看出,本文給出的全自由度模型分析結果與傳統有限元結果較為吻合.而解析解[2]由低階傅里葉展開得到,溫度按照三角函數沿圓周分布,故而與數值解之間存在一定誤差.圖6 給出了FOM-TANCF 和ROM-TANCF 的溫度分布.圖7 給出了FOM-TANCF 和ROM-TANCF 之間的溫度誤差.溫度均方根誤差(TRMS),用來評估所提出的縮減后的模型和全自由度模型之間溫度的誤差,表達方式如下

圖5 不同時刻圓周溫度分布對比Fig.5 Comparison of section temperature distribution at different times

其中n是節點個數.

溫度坐標的總自由度為30,ROM-TANCF 中保留了10 個溫度坐標.時間步長為0.01 s.計算時間從321.3 s 減少到119.8 s,降幅為62.71%.如圖5所示,FOM-TANCF 和FEM 的結果一致.因此,FOM-TANCF 的結果作為基準.如圖6 和圖7 所示,FOM-TANCF 和ROM-TANCF 的計算結果吻合良好,證明該方法在保證足夠精度的前提下能有效提高求解傳熱方程的效率.

圖6 FOM-TANCF 和ROM-TANCF 溫度分布對比Fig.6 Comparison of section temperature distribution between FOMTANCF and ROM-TANCF

圖7 FOM-TANCF 和ROM-TANCF 溫度誤差Fig.7 Error of temperature between FOM-TANCF and ROM-TANCF

5.2 薄板熱膨脹

本節提出了一個長方形薄板模型來測試熱膨脹.如圖8(a)所示.將系統在長度方向上劃分為6 個熱柔耦合薄板單元.系統分為兩個子結構,每個子結構劃分為3 個單元,如圖8(b)所示.初始溫度為0 K.總仿真時間為6 s.其中前5 s 是恒定的內熱,最后1 s 內熱停止,其他參數見表2.均方根誤差(RMS),用來評估所提出的縮減后的模型和全自由度模型之間的誤差,表達方式如下

圖8 薄板熱膨脹模型Fig.8 Thermal deformation of thin plate mode

表2 薄板熱膨脹板模型參數Table 2 Parameters of the thermal expansion thin plate

給定的閾值 δk和 δT由式(24)獲得,分別為2.5 mm和0.01 K.時間步長為0.1 ms.計算時間從5089.7 s減少到1422.5 s,減少了72.05%.動力學自由度從126 個減少到70 個,其中子結構1 保留26 個,子結構2 保留26 個,子結構連接邊界約束分別為18 個.為了方便起見,可以將記為70 (26+26+18).同樣,對于傳熱方程,自由度從42 個減少到26 個(10+10 +6).采用ANSYS 中的SOLID90單元(網格劃分為200×50×2),采用間接耦合求解.如圖9 和10 所示,FOM-TANCF,ROM-TANCF 和ANSYS 的結果一致.圖11 中給出了RMS 誤差.計算結果表明,該方法能在有效提高求解熱膨脹問題的計算效率的同時保證較高的計算精度.

圖9 檢測點溫度Fig.9 The test point temperature

圖10 檢測點在長度方向位置Fig.10 The test point position in length direction

圖11 FOM-TANCF 和ROM-TANCF 誤差Fig.11 Error of deformation between FOM-TANCF and ROM-TANCF

5.3 柔性太陽能電池板

如圖12(a)所示,最初水平放置的柔性太陽能電池板一端鉸接,并在微重力的作用下下落.太陽能電池板由中間的薄板結構和兩側的梁結構組成.薄板單元和梁單元節點通過約束方程連接.太陽熱流水平射入,初始時不受熱流照射,隨著微重力下擺時,受熱面積逐漸增大,溫度逐漸升高.薄板結構的網格為4×2 (長×寬).每個梁離散為8 個圓形截面ANCF梁單元[39].系統分為兩個子結構,如圖12(b)所示.總模擬時間設置為10 s.初始溫度為0 K,柔性太陽能電池板的其他參數如表3 所示.

圖12 柔性太陽能電池板模型Fig.12 The flexible solar panel model

表3 柔性太陽能電池板模型參數Table 3 Parameters of the flexible solar panel

給定的閾值 δk和δT分別為1 mm 和0.01 K.時間步長為0.2 ms.計算時間從30 363.6 s 減少到9293.2 s,減少了69.39%.動力學的自由度從351 個減少到227 個(100+100+27).對于傳熱方程,自由度從117 個減少到75 個(33+33+9).

X方向上的測試點位置如圖13 所示.梁和薄板的熱膨脹系數不同,如表3 所示.由于大變形和旋轉,柔性太陽能電池板的下擺導致不均勻溫升和熱應變.與沒有溫度變化的情況相比,不均勻的熱應變會導致測試點軌跡的不同.RMS 誤差如圖14 所示,FOM-TANCF 和ROM-TANCF 的計算結果一致.圖15給出了由ROM-TANCF 獲得的柔性太陽能電池板在不同時間的運動過程和溫度分布.計算結果表明,該方法能有效提高大變形、大轉動的熱柔耦合系統的效率.

圖13 檢測點X-方向位置對比Fig.13 Comparison of test point position in X-direction

圖14 FOM-TANCF 和ROM-TANCF 誤差Fig.14 Error of deformation between FOM-TANCF and ROM-TANCF

圖15 柔性太陽能電池板不同時刻構型Fig.15 Motion process of the flexible solar panel

5.4 剛柔熱耦合拋物線天線

本節使用剛柔熱耦合天線驗證所提出算法的有效性.剛柔熱耦合天線模型如圖16 所示.采用薄板單元模擬拋物面天線,三維梁單元模擬拋物面天線上的肋梁.兩種單元在公共節點處通過約束方程連接.如圖16 所示,柔性天線分為兩個子結構.ANCF參考節點用于建模中心剛性輪轂.中心輪轂的質量為300 kg,中心輪轂繞主軸的轉動慣量分別為Jxx=100 kg·m2,Jyy=100 kg·m2,Jzz=500 kg·m2.太陽熱流的入射方向為[5,0,-2].仿真時長為50 s,初始溫度為0 K,其他參數見表4.驅動扭矩施加在中心剛性輪轂上,驅動力矩表達式如下

表4 剛柔熱耦合天線模型參數Table 4 Parameters of the rigid-flexible-thermal coupling satellite

圖16 剛柔熱耦合天線模型Fig.16 The rigid-flexible-thermal coupling satellite schematic model

閾值 δk和 δT分別為3 mm 和0.01 K.時間步長為0.1 ms.測試點A和B在Y方向上的位置分量如圖17 所示,溫度如圖18 所示.RMS 誤差如圖19 所示.不同時間的溫度分布如圖20 所示.可以看出相較于直徑達10 m 的天線,RMS 最大值不超過3 cm,與此同時,計算時間從337 303.2 s 減少到67 466.5 s,減少了79.99%.動力學的自由度從708 個減少到554 個(250+250+54).對于傳熱方程,自由度從232 個減少到148 個(65+65+18).盡管本文提出的方法需要在位移或溫度增量過大時,對子結構進行重新線性化,更新模態基,但是ANCF全自由度模型每個時間步長內彈性力及其雅可比求解都非常耗時,降階模型從整體上體現出更高的仿真求解效率.

圖17 A 和B 點Y 方向位置對比Fig.17 Comparison of test points A and B position in Y-direction

圖18 A 和B 點溫度Fig.18 Temperature of test points A and B

圖19 FOM-TANCF 和 ROM-TANCF 誤差Fig.19 Error of deformation between FOM-TANCF and ROM-TANCF

圖20 不同時刻溫度分布Fig.20 Temperature distribution at different times

6 結論

本文提出了一種基于改進模態綜合法的剛柔熱耦合多體系統模型降階方法.ANCF 單元形函數被同時用于描述具有大變形和大范圍運動的多柔體系統及柔性體內溫度場差值離散.采用ANCF 參考節點用于描述剛體運動,得到了統一描述的剛柔熱耦合多體系統模型.采用泰勒展開法對系統動力學和傳熱方程進行線性化處理,使得線性化區間內,柔性體的切線剛度矩陣可視為常值矩陣.采用了改進的模態綜合法對線性化后的動力學方程和傳熱方程進行降階.子結構之間通過約束方程連接,以確保連接精度及邊界連續性.最后,通過四個數值算例表明,降階后的模型計算結果與全自由度模型計算結果一致,并可以有效提高計算效率.半圓形薄板的純熱傳導算例考慮了熱流輸入和表面自熱輻射.計算結果與以往文獻中的結果吻合,計算時間減少了62.71%,保留了33.33%的自由度.薄板的熱膨脹算例中,計算結果與ANSYS 計算結果吻合較好,驗證了該方法的準確性.對于柔性太陽能電池板和剛柔熱耦合天線,計算時間可分別減少69.39%和79.99%.算例結果表明,該方法能有效地提高剛柔熱耦合系統的計算效率,并保持較高的精度.此外,該模型降階方法可以推廣到阻尼系統和其他多個物理耦合的領域.

附錄

與熱應變對應的廣義力可以寫成

其中,ε 是應變向量,Eε是材料彈性系數矩陣,εH是熱應變向量.熱應變對應的廣義力的雅可比矩陣可以寫為

式中,α 為熱膨脹系數,ΔT為溫度增量,薄板單元材料彈性系數矩陣Eε可以在文獻[18]中找到.

式中,α 為熱膨脹系數,ΔT為溫度增量,三維梁單元材料彈性系數矩陣Eε可以在文獻[34]中找到.