熱-動(dòng)力學(xué)耦合多體系統(tǒng)建模與降階1)

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

* (哈爾濱工業(yè)大學(xué)機(jī)電工程學(xué)院,哈爾濱 150001)

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

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

?? (西安建筑科技大學(xué)西部綠色建筑國(guó)家重點(diǎn)實(shí)驗(yàn)室,西安 710055)

引言

太陽(yáng)能電池板和拋物面天線等典型的航天器板結(jié)構(gòu)具有重量輕、尺寸大的特點(diǎn).此類航天器的動(dòng)力學(xué)分析與控制問題越來越多地需要考慮到空間環(huán)境溫度的影響.由太陽(yáng)輻射導(dǎo)致的極端熱載荷激勵(lì)引起位移場(chǎng)與溫度場(chǎng)的強(qiáng)耦合響應(yīng),進(jìn)而導(dǎo)致熱致振動(dòng)或熱跳變[1-4].因此,剛?cè)釤狁詈舷到y(tǒng)的精確建模對(duì)航天器的安全至關(guān)重要.

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

隨著ANCF 計(jì)算規(guī)模增大,數(shù)值效率低的問題日漸突出,模型降階成為了重要的研究方向.在文獻(xiàn)[19-20]的研究中,本征正交分解(POD)方法被用于ANCF 單元的模型降階.然而,降階模型必須通過實(shí)驗(yàn)數(shù)據(jù)庫(kù)或計(jì)算全階模型仿真結(jié)果來獲得.模態(tài)綜合法(CMS) 可以避免這一缺點(diǎn).Kobayashi 等[21]在軸向變形較小的前提下,使用固定界面模態(tài)綜合法(也稱Craig-Bampton 方法)來縮減ANCF 平面梁?jiǎn)卧５亩囿w系統(tǒng).Sun 等[22-23]使用自由界面模態(tài)綜合法來提高帶滑動(dòng)鉸的柔性梁的計(jì)算效率.該方法同樣采用了相同的軸向小變形的假設(shè).Tang 等[24]將整個(gè)柔性體運(yùn)動(dòng)過程劃分為若干子區(qū)間并在每個(gè)子區(qū)間內(nèi)對(duì)系統(tǒng)方程進(jìn)行線性化處理,從而使得切線剛度陣定常,基于模態(tài)截?cái)嗟慕惦A方法得以應(yīng)用.此外作者也采用了Craig-Bampton 方法對(duì)ANCF 梁?jiǎn)卧枋龅亩嗳狍w系統(tǒng)進(jìn)行子結(jié)構(gòu)劃分,稱之為CB-ANCF 方法.在此基礎(chǔ)上,作者將此方法推廣到了具有阻尼的多柔體系統(tǒng),以減少黏彈性多柔體系統(tǒng)的自由度[25].基于文獻(xiàn)[24-25] 的研究,Tian等[26]使用自由界面模態(tài)綜合法來減少ANCF 梁?jiǎn)卧枋龅亩嗳狍w系統(tǒng)的自由度,即F-ANCF 方法.

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

本文提出了一種剛?cè)釤狁詈舷到y(tǒng)模型降階的方法.引入溫度梯度,采用高階插值對(duì)溫度場(chǎng)進(jìn)行離散,使位移場(chǎng)和溫度場(chǎng)采用統(tǒng)一的ANCF 單元形函數(shù)進(jìn)行插值離散,避免數(shù)值求解中網(wǎng)格錯(cuò)配和誤差映射的問題.ANCF 參考節(jié)點(diǎn)用于描述中心剛體.采用泰勒展開對(duì)動(dòng)力學(xué)方程和傳熱方程進(jìn)行線性化處理,得到常值切線剛度陣.引入改進(jìn)的模態(tài)綜合法分別對(duì)動(dòng)力學(xué)方程和傳熱方程進(jìn)行自由度縮減.最后,以剛?cè)釤狁詈洗罂趶綊佄锩嫣炀€為例驗(yàn)證了該方法的有效性.

1 熱-動(dòng)力學(xué)耦合建模

本文擬分別采用熱耦合的ANCF 薄板[18]與三維梁?jiǎn)卧猍17]分別構(gòu)造星載拋物面天線及其肋梁模型.而衛(wèi)星本體視為剛體,由ANCF 參考節(jié)點(diǎn)[32]描述.由此可將衛(wèi)星剛?cè)釤狁詈舷到y(tǒng)納入ANCF 框架內(nèi)描述和建模.由于采用了全局坐標(biāo)系下定義的位置與梯度向量,單元內(nèi)部任意點(diǎn)的空間位置可以定義為r=Se,S為單元形函數(shù)矩陣,e為單元節(jié)點(diǎn)向量.ANCF 薄板單元基于Kirchhoff-Love 板理論,其節(jié)點(diǎn)向量包含中面四個(gè)角點(diǎn)的位置和平面內(nèi)梯度向量,即其中上標(biāo)P 代指板單元,i=1,2,3,4 為節(jié)點(diǎn)號(hào).ANCF 三維梁?jiǎn)卧?jié)點(diǎn)為中軸線兩端點(diǎn),每個(gè)節(jié)點(diǎn)使用全局位置向量和沿三個(gè)方向的梯度向量作為節(jié)點(diǎn)坐標(biāo),可以表示為:eB=上標(biāo)B 代指梁?jiǎn)卧?j=1,2.

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

其中,M為系統(tǒng)質(zhì)量矩陣,Qe是彈性力.薄板單元和梁?jiǎn)卧膹椥粤Ψ謩e在文獻(xiàn)[33-34]中給出.QH是對(duì)應(yīng)于熱應(yīng)變的廣義力,表達(dá)式在附錄中給出.Ce為約束方程,?Ce/?e則為約束方程對(duì)廣義坐標(biāo)e的偏導(dǎo)數(shù),λ 是拉格朗日乘子,Qext代表廣義外力.

由于連續(xù)體內(nèi)溫度場(chǎng)采用了與運(yùn)動(dòng)學(xué)描述相同的離散格式,單元節(jié)點(diǎn)處的溫度及其梯度被用作傳熱方程的廣義坐標(biāo),如圖1 所示.薄板單元的溫度廣義坐標(biāo)為其中Ti(i=1,2,3,4)為節(jié)點(diǎn)i的溫度,Txi和Tyi分別表示節(jié)點(diǎn)i在x和y方向上的溫度梯度.同理,三維梁?jiǎn)卧臏囟葟V義坐標(biāo)可以寫成

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

與位置插值同理,單元內(nèi)任意點(diǎn)的溫度可以寫為

其中,ST是溫度場(chǎng)的單元形狀函數(shù),其元素與描述位移場(chǎng)形狀函數(shù)相同,eT表示單元的節(jié)點(diǎn)廣義溫度坐標(biāo)向量.同樣地,上角標(biāo)P 和B 分別指代薄板單元和梁?jiǎn)卧?

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

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

2 動(dòng)力學(xué)和傳熱方程的線性化

基于ANCF 建立的柔性體動(dòng)力學(xué)方程具有常值質(zhì)量陣的優(yōu)勢(shì),但其彈性力卻是節(jié)點(diǎn)坐標(biāo)向量的高度非線性函數(shù),切線剛度陣亦然.這就導(dǎo)致基于模態(tài)截?cái)嗟慕惦A方法無法直接應(yīng)用.為此,將系統(tǒng)的整個(gè)運(yùn)動(dòng)過程劃分為多個(gè)區(qū)間,認(rèn)為柔性體在每個(gè)區(qū)間內(nèi)部的動(dòng)力學(xué)行為是準(zhǔn)靜態(tài)的.利用一階泰勒展開將基于初始構(gòu)型定義的全量平衡方程變換為基于當(dāng)前區(qū)間起始構(gòu)型的增量平衡形式.每個(gè)區(qū)間內(nèi)的切線剛度矩陣可視為定常,進(jìn)而采用改進(jìn)的模態(tài)綜合法對(duì)系統(tǒng)進(jìn)行降階.

根據(jù)泰勒展開式,每個(gè)準(zhǔn)靜態(tài)區(qū)間的初始構(gòu)型處的動(dòng)力學(xué)方程可展開為

其中,e0是每個(gè)準(zhǔn)靜態(tài)區(qū)間初始構(gòu)型處的廣義坐標(biāo).同理,eT0是初始廣義溫度坐標(biāo).JQe是彈性力雅可比矩陣(表達(dá)式見文獻(xiàn)[33-34]).JQH為熱應(yīng)變廣義力的雅可比矩陣,其表達(dá)式在附錄中給出.約束方程也可以線性化為

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

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

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

傳熱約束方程可展開為

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

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

3 改進(jìn)的模態(tài)綜合法

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

圖2 系統(tǒng)子結(jié)構(gòu)劃分Fig.2 The mesh of the system and substructures

在動(dòng)力學(xué)方程部分,子結(jié)構(gòu)j的動(dòng)力學(xué)方程齊次形式如下

其中,[]j表示子結(jié)構(gòu)j.通過模態(tài)分析獲得當(dāng)前子結(jié)構(gòu)的特征值和特征向量

其中,ωn是第n個(gè)固有頻率,[Φ]j為子結(jié)構(gòu)振型,是子結(jié)構(gòu)j的自由度數(shù).

子結(jié)構(gòu)j的廣義坐標(biāo)增量可以表示為

以兩個(gè)子結(jié)構(gòu)組成的系統(tǒng)為例,其廣義坐標(biāo)可以表示為

顯然,每個(gè)子結(jié)構(gòu)的坐標(biāo)是獨(dú)立的,除了連接處的自由度.在本研究中,直接施加界面位置和梯度的約束方程,以確保子結(jié)構(gòu)連接處的精度,即

將式(14)和式(15)代入式(6),系統(tǒng)動(dòng)力學(xué)方程可以表示為

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

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

式中,ωm為第m熱模態(tài)頻率,物理意義為熱衰減系數(shù)[27],[Ψ]j是對(duì)應(yīng)于熱特征值的特征向量是子結(jié)構(gòu)j的溫度坐標(biāo)數(shù).傳熱方程的特征值的求解在文獻(xiàn)[36]中給出.

子結(jié)構(gòu)j的廣義溫度坐標(biāo)增量可以表示為

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

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

4 求解策略

為準(zhǔn)確模擬熱-動(dòng)力學(xué)耦合效應(yīng)對(duì)系統(tǒng)的影響,需要同步求解傳熱方程和動(dòng)力學(xué)方程.本文采用廣義-α 積分來求解傳熱方程和動(dòng)力學(xué)方程.圖3 給出了算法的流程圖.算法的數(shù)值參數(shù) αf,αm,β ,γ 等的取值見參考文獻(xiàn)[37-38].

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

在第2 節(jié)中,使用泰勒展開法將動(dòng)力學(xué)方程和傳熱方程線性化,將非線性方程轉(zhuǎn)化為線性方程.需要指出的是,線性方程等式成立的前提是廣義位置坐標(biāo)的增量應(yīng)該很小.當(dāng)?shù)趈個(gè)子結(jié)構(gòu)的最大位移增量過大時(shí),第j子結(jié)構(gòu)應(yīng)重新線性化.同樣,溫度的升高也會(huì)引起熱應(yīng)變的變化,同樣會(huì)影響切線剛度陣,因此,當(dāng)?shù)趈子結(jié)構(gòu)的溫度增量過大時(shí),第j子結(jié)構(gòu)也應(yīng)重新線性化.保證 δk和 δT足夠小,將累計(jì)誤差控制在一定范圍內(nèi).本文給出閾值的經(jīng)驗(yàn)計(jì)算公式

其中L是梁結(jié)構(gòu)的長(zhǎng)度或者薄板結(jié)構(gòu)長(zhǎng)度.通常,對(duì)于動(dòng)力學(xué)方程,保留的自由度nred建議如下nred=ntal/3～4ntal/5.對(duì)于傳熱方程,保留的自由度建議如下

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

5 數(shù)值算例

本節(jié)給出了四個(gè)數(shù)值算例來驗(yàn)證所提出方法的準(zhǔn)確性和有效性.為方便起見,完整自由度模型記為FOM-TANCF,降階模型記為ROM-TANCF.

5.1 半圓薄板的純導(dǎo)熱

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

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

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

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

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

其中n是節(jié)點(diǎn)個(gè)數(shù).

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

圖6 FOM-TANCF 和ROM-TANCF 溫度分布對(duì)比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 薄板熱膨脹

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

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

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

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

圖9 檢測(cè)點(diǎn)溫度Fig.9 The test point temperature

圖10 檢測(cè)點(diǎn)在長(zhǎng)度方向位置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 柔性太陽(yáng)能電池板

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

圖12 柔性太陽(yáng)能電池板模型Fig.12 The flexible solar panel model

表3 柔性太陽(yáng)能電池板模型參數(shù)Table 3 Parameters of the flexible solar panel

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

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

圖13 檢測(cè)點(diǎn)X-方向位置對(duì)比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 柔性太陽(yáng)能電池板不同時(shí)刻構(gòu)型Fig.15 Motion process of the flexible solar panel

5.4 剛?cè)釤狁詈蠏佄锞€天線

本節(jié)使用剛?cè)釤狁詈咸炀€驗(yàn)證所提出算法的有效性.剛?cè)釤狁詈咸炀€模型如圖16 所示.采用薄板單元模擬拋物面天線,三維梁?jiǎn)卧M拋物面天線上的肋梁.兩種單元在公共節(jié)點(diǎn)處通過約束方程連接.如圖16 所示,柔性天線分為兩個(gè)子結(jié)構(gòu).ANCF參考節(jié)點(diǎn)用于建模中心剛性輪轂.中心輪轂的質(zhì)量為300 kg,中心輪轂繞主軸的轉(zhuǎn)動(dòng)慣量分別為Jxx=100 kg·m2,Jyy=100 kg·m2,Jzz=500 kg·m2.太陽(yáng)熱流的入射方向?yàn)閇5,0,-2].仿真時(shí)長(zhǎng)為50 s,初始溫度為0 K,其他參數(shù)見表4.驅(qū)動(dòng)扭矩施加在中心剛性輪轂上,驅(qū)動(dòng)力矩表達(dá)式如下

表4 剛?cè)釤狁詈咸炀€模型參數(shù)Table 4 Parameters of the rigid-flexible-thermal coupling satellite

圖16 剛?cè)釤狁詈咸炀€模型Fig.16 The rigid-flexible-thermal coupling satellite schematic model

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

圖17 A 和B 點(diǎn)Y 方向位置對(duì)比Fig.17 Comparison of test points A and B position in Y-direction

圖18 A 和B 點(diǎn)溫度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 不同時(shí)刻溫度分布Fig.20 Temperature distribution at different times

6 結(jié)論

本文提出了一種基于改進(jìn)模態(tài)綜合法的剛?cè)釤狁詈隙囿w系統(tǒng)模型降階方法.ANCF 單元形函數(shù)被同時(shí)用于描述具有大變形和大范圍運(yùn)動(dòng)的多柔體系統(tǒng)及柔性體內(nèi)溫度場(chǎng)差值離散.采用ANCF 參考節(jié)點(diǎn)用于描述剛體運(yùn)動(dòng),得到了統(tǒng)一描述的剛?cè)釤狁詈隙囿w系統(tǒng)模型.采用泰勒展開法對(duì)系統(tǒng)動(dòng)力學(xué)和傳熱方程進(jìn)行線性化處理,使得線性化區(qū)間內(nèi),柔性體的切線剛度矩陣可視為常值矩陣.采用了改進(jìn)的模態(tài)綜合法對(duì)線性化后的動(dòng)力學(xué)方程和傳熱方程進(jìn)行降階.子結(jié)構(gòu)之間通過約束方程連接,以確保連接精度及邊界連續(xù)性.最后,通過四個(gè)數(shù)值算例表明,降階后的模型計(jì)算結(jié)果與全自由度模型計(jì)算結(jié)果一致,并可以有效提高計(jì)算效率.半圓形薄板的純熱傳導(dǎo)算例考慮了熱流輸入和表面自熱輻射.計(jì)算結(jié)果與以往文獻(xiàn)中的結(jié)果吻合,計(jì)算時(shí)間減少了62.71%,保留了33.33%的自由度.薄板的熱膨脹算例中,計(jì)算結(jié)果與ANSYS 計(jì)算結(jié)果吻合較好,驗(yàn)證了該方法的準(zhǔn)確性.對(duì)于柔性太陽(yáng)能電池板和剛?cè)釤狁詈咸炀€,計(jì)算時(shí)間可分別減少69.39%和79.99%.算例結(jié)果表明,該方法能有效地提高剛?cè)釤狁詈舷到y(tǒng)的計(jì)算效率,并保持較高的精度.此外,該模型降階方法可以推廣到阻尼系統(tǒng)和其他多個(gè)物理耦合的領(lǐng)域.

附錄

與熱應(yīng)變對(duì)應(yīng)的廣義力可以寫成

其中,ε 是應(yīng)變向量,Eε是材料彈性系數(shù)矩陣,εH是熱應(yīng)變向量.熱應(yīng)變對(duì)應(yīng)的廣義力的雅可比矩陣可以寫為

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

式中,α 為熱膨脹系數(shù),ΔT為溫度增量,三維梁?jiǎn)卧牧蠌椥韵禂?shù)矩陣Eε可以在文獻(xiàn)[34]中找到.