張佳明，王文瑞，馬新杰，溫曉東

(1.北京科技大學(xué) 機(jī)械工程學(xué)院， 北京 100083; 2.流體與材料相互作用教育部重點實驗室， 北京 100083;3.天津航天機(jī)電設(shè)備研究所 天津市宇航智能裝備技術(shù)企業(yè)重點實驗室， 天津 300301)

近空間飛行器擁有快速響應(yīng)、突防能力強(qiáng)等特點，能夠在短時間內(nèi)精準(zhǔn)打擊長距離軍事目標(biāo)，在國防科研、空間技術(shù)研究中具有極其重要的戰(zhàn)略意義[1-2]。飛行器在高馬赫數(shù)飛行過程中，周圍空氣受到激波的強(qiáng)烈壓縮，對飛行器表面產(chǎn)生劇烈的氣動熱力載荷。舵翼結(jié)構(gòu)作為飛行器關(guān)鍵部件，為了確保燃料供給能夠滿足高速、長距離飛行，要求盡量減小結(jié)構(gòu)尺寸，結(jié)構(gòu)剛度的降低導(dǎo)致更大的變形量，影響其氣動外形，進(jìn)一步增強(qiáng)氣動熱、氣動力與結(jié)構(gòu)響應(yīng)之間的耦合效應(yīng)[3-5]，直接影響飛行穩(wěn)定性，嚴(yán)重時造成偏離飛行軌道甚至掉落，影響精確打擊效果，這也是目前近空間飛行器熱防護(hù)設(shè)計與結(jié)構(gòu)設(shè)計過程中必然要面對的問題。

為此，多位研究者開展了相關(guān)的試驗項目。夏吝時等[6]對主動式熱疏導(dǎo)簡化前緣金屬開展地面熱環(huán)境模擬試驗。吳大方等[7-8]對近空間飛行器復(fù)合材料先后開展了1 100 ℃和1 500 ℃的熱振聯(lián)合實驗。蔡禮港、歐朝、李強(qiáng)等[9-11]在電弧風(fēng)洞、激波風(fēng)洞中對飛行器結(jié)構(gòu)與熱防護(hù)材料開展了實驗研究。然而，開展試驗周期長、難度大、費用高，使用數(shù)值仿真方法對近空間飛行器結(jié)構(gòu)溫度、變形進(jìn)行預(yù)測，能夠更為高效地為熱防護(hù)設(shè)計與結(jié)構(gòu)設(shè)計提供指導(dǎo)。李佳偉等[12]采用一體化求解方法，將流場與結(jié)構(gòu)溫度場作為一個物理場，采用統(tǒng)一的控制方程進(jìn)行求解，避開了傳統(tǒng)耦合方法在迭代計算時場之間的大量數(shù)據(jù)交換。季衛(wèi)棟等[13]采用統(tǒng)一的積分方程組作為氣動加熱和結(jié)構(gòu)熱物理過程的控制方程，對物理場進(jìn)行統(tǒng)一的迎風(fēng)格式有限體積方法離散，對二維圓管模型的氣動加熱和結(jié)構(gòu)傳熱問題進(jìn)行了模擬。陳梅潔等[14]通過順序耦合方法，基于Transitionk-k1-ω湍流模型計算得到了溫度場和應(yīng)力場的分布。梁杰等[15]采用蒙特卡洛方法，通過分區(qū)求解策略對大型航天器離軌再入隕落過程中的氣動力和氣動熱特性進(jìn)行了數(shù)值模擬。

當(dāng)前研究者更多關(guān)注于多物理場耦合方法的實現(xiàn)，而在實現(xiàn)仿真的基礎(chǔ)上，開展各物理場之間的耦合關(guān)系分析，對于熱防護(hù)材料設(shè)計及其結(jié)構(gòu)設(shè)計有著更為現(xiàn)實的意義[16]。為此，本文通過分區(qū)耦合方法，開展舵翼結(jié)構(gòu)流-熱-固耦合仿真，對不同飛行迎角狀態(tài)下結(jié)構(gòu)變形場、溫度場進(jìn)行計算，實現(xiàn)氣動力、氣動熱及結(jié)構(gòu)耦合響應(yīng)的量化分析。

1 多場耦合理論

1.1 數(shù)值計算理論

流-熱-固多場耦合根據(jù)以下6個方程進(jìn)行數(shù)值計算。

1) 質(zhì)量守恒方程。

(1)

式(1)中：▽代表拉普拉斯算子；V為速度矢量(u、v、w)，指流體微元在x、y、z等3個方向的速度分量(m/s)；ρf為流體密度(kg/m3)；t為時間(s)。

2) 動量守恒方程。

(2)

3) 能量守恒方程。

(3)

式(3)中：E為總能量(J)；ρ為材料密度(kg/m3)；fw=fxu+fyv+fzw為單位質(zhì)量的流體體積力所做的功(J)。

4) 氣體狀態(tài)方程。

p=Z(R*/M)ρfT

(4)

式(4)中：M表示氣體摩爾質(zhì)量(kg/mol)；R*為摩爾氣體常數(shù)，R*=8.314 5 J/(mol·K)；T為溫度場(K)；Z為壓縮因子。

5) 傳熱方程。

(5)

式(5)中：λ為材料的導(dǎo)熱系數(shù)(W/(m·K))；C為比熱容(J/kg·K)；Q為物體內(nèi)部的熱源強(qiáng)度(W/kg)。

6) 結(jié)構(gòu)響應(yīng)方程。

[σ]=[D][ε]

(6)

式(6)中，[σ]為應(yīng)力矩陣；[D]為結(jié)構(gòu)材料彈性矩陣；[ε]為結(jié)構(gòu)的總應(yīng)變矩陣。

1.2 多物理場耦合方法

舵翼結(jié)構(gòu)在近空間環(huán)境中飛行時，氣動加熱過程是一個持續(xù)的、非穩(wěn)態(tài)的過程，并且是流場、溫度場和結(jié)構(gòu)場相互耦合作用的過程，因此采用流-熱-固耦合方法進(jìn)行多物理場耦合問題的求解。熱流固耦合模型如圖1所示，在流體域內(nèi)使用基于有限體積法的Fluent求解器，對流場的熱流和壓力進(jìn)行求解，使用System Coupling模塊，通過流固耦合界面將熱載與氣動力載荷傳遞至固體域。在固體域內(nèi)使用Transient Structural求解器，對固體域的溫度與變形進(jìn)行求解，使用System Coupling模塊通過流固耦合界面將溫度場與變形場傳遞至流體域，繼續(xù)進(jìn)行流體域的下一步計算，直到達(dá)到設(shè)置的耦合計算時間。

圖1 熱流固耦合模型示意圖

2 仿真算法校驗

為驗證本文選用多物理場耦合算法的準(zhǔn)確性，利用經(jīng)典的圓管繞流試驗作為驗證算例。該試驗是1987年美國科學(xué)家 Allan R.Wieting在NASA LANGLEY 8-foot高溫風(fēng)洞中進(jìn)行的，試驗如圖2所示。

圖2 圓管繞流風(fēng)洞試驗示意圖

風(fēng)洞試驗中，不銹鋼圓管是完全對稱的圓柱體，因此本文中采用1/2模型進(jìn)行計算，以提高計算效率。利用SolidWorks軟件建立流場及圓柱模型，采用結(jié)構(gòu)化網(wǎng)格，在邊界層處和激波位置處進(jìn)行網(wǎng)格加密處理。流體域和固體域在耦合面上網(wǎng)格節(jié)點是連續(xù)的，保證流固耦合面上的幾何相容性，模型網(wǎng)格如圖3所示。流體域網(wǎng)格數(shù)量為22 866，固體域網(wǎng)格數(shù)量為3 249，網(wǎng)格質(zhì)量均大于0.985。流體域入口、出口設(shè)置為遠(yuǎn)場邊界條件，流場參數(shù)設(shè)置如表1所示，耦合面按流固耦合邊界條件設(shè)置。圓管內(nèi)徑為25.4 mm，外徑為38.1 mm，材料為1Cr18Ni9Ti不銹鋼，材料參數(shù)如表2所示，圓管初始溫度為294.4 K。

圖3 網(wǎng)格邊界條件示意圖

表1 來流參數(shù)

表2 不銹鋼材料參數(shù)

仿真耦合求解時間設(shè)置為2s，與文獻(xiàn)[17]相同。圖4為仿真流場溫度分布云圖與文獻(xiàn)試驗激波紋影，二者形狀相吻合。圖5表示了仿真結(jié)果與文獻(xiàn)中的流場溫度突變位置、流場溫度，二者分別相差0.071%與2.369%。圖6表示了仿真結(jié)果與文獻(xiàn)中的圓管結(jié)構(gòu)溫度分布、應(yīng)力分布，二者分別相差0.345%和0.567%，具體數(shù)據(jù)如表3所示。由上述分析可知，仿真結(jié)果與文獻(xiàn)所述基本吻合，說明本文所使用的多場耦合方法是準(zhǔn)確有效的。

圖4 外流場參數(shù)分布圖

圖5 駐點徑向流場參數(shù)曲線

圖6 圓管結(jié)構(gòu)溫度與應(yīng)力分布圖

表3 仿真與文獻(xiàn)數(shù)據(jù)對比

3 舵翼結(jié)構(gòu)模型耦合計算結(jié)果與分析

3.1 有限元建模

舵翼結(jié)構(gòu)尺寸如圖7(a)所示，材料為D6AC高溫合金鋼，材料參數(shù)如表4所示，計算模型流體域、固體域網(wǎng)格分別如圖7(b)、圖7(c)所示。模型計算域的邊界條件如圖7(d)所示，壁面邊界選擇無滑移邊界條件，將初始壁溫設(shè)置為216.65 K，流體域外邊界設(shè)置為壓力遠(yuǎn)場邊界條件，模擬飛行高度為20 km，來流速度5Ma，壓強(qiáng)5 529.3 Pa，溫度216.65 K，密度0.088 91 kg/m3。固體結(jié)構(gòu)場設(shè)置為彈性結(jié)構(gòu)，尾部施加固定約束，對稱面添加對稱約束，將舵翼結(jié)構(gòu)與外流場交界面設(shè)置為流-固耦合面，在界面上進(jìn)行分區(qū)耦合之間的數(shù)據(jù)傳輸。根據(jù)舵翼結(jié)構(gòu)飛行特點，設(shè)置飛行迎角為-10°、-5°、0°、5°、10°進(jìn)行計算，耦合時間為40 s。

圖7 舵翼結(jié)構(gòu)與計算域網(wǎng)格

表4 D6AC鋼材料參數(shù)

為了分析舵翼結(jié)構(gòu)在近空間飛行狀態(tài)下的氣動力、氣動熱與結(jié)構(gòu)響應(yīng)之間的關(guān)系，將仿真工況設(shè)置為耦合計算(考慮氣動力、氣動力耦合效應(yīng))、僅氣動熱計算(僅考慮熱變形)、僅氣動力計算(僅考慮氣動力變形)和無結(jié)構(gòu)變形因素計算4種工況。

3.2 結(jié)構(gòu)變形分析

舵翼結(jié)構(gòu)在X、Y、Z3個方向發(fā)生變形，在不同迎角狀態(tài)下發(fā)生拉壓、彎曲和熱膨脹變形，結(jié)構(gòu)在Z方向只發(fā)生熱膨脹變形而無拉壓和彎曲變形，因此只對舵翼結(jié)構(gòu)X方向(縱向)和Y方向(橫向)的變形進(jìn)行分析。

對舵翼結(jié)構(gòu)施加耦合工況，圖8為不同迎角狀態(tài)下的舵翼結(jié)構(gòu)縱向變形云圖，圖9為舵翼結(jié)構(gòu)不同位置的縱向變形曲線，其中橫坐標(biāo)為舵翼結(jié)構(gòu)所處位置，X=0處為前緣頂端位置，如圖7(a)所示。

圖8 舵翼結(jié)構(gòu)X方向變形云圖

圖9 舵翼結(jié)構(gòu)不同迎角下縱向變形曲線

舵翼結(jié)構(gòu)在近空間飛行狀態(tài)下發(fā)生氣動加熱現(xiàn)象，使舵翼結(jié)構(gòu)產(chǎn)生縱向變形。由于舵翼結(jié)構(gòu)尾部采用固定約束，類似于懸臂梁受載方式，當(dāng)迎角為-10°時，舵翼結(jié)構(gòu)尖角的平分線與來流方向夾角最小，此時舵翼結(jié)構(gòu)彎曲撓度最小，而縱向變形量最大。隨著迎角增大撓度隨之增加，縱向變形量逐漸減小。

圖10為舵翼結(jié)構(gòu)不同位置的橫向變形曲線，與縱向變形規(guī)律不同，隨著迎角增大，橫向變形量逐漸增大。舵翼結(jié)構(gòu)尾部采用固定約束，受到橫向氣動力作用，發(fā)生彎曲變形。另一方面，由于氣動加熱使結(jié)構(gòu)發(fā)生熱膨脹，迎風(fēng)側(cè)溫度較高，膨脹量大，背風(fēng)側(cè)溫度低，膨脹量小，造成舵翼結(jié)構(gòu)沿背風(fēng)側(cè)發(fā)生彎曲。在氣動力和氣動熱的共同作用下，舵翼結(jié)構(gòu)發(fā)生彎曲變形，迎角越大，上述現(xiàn)象越明顯，橫向變形量也就越大。

圖10 舵翼結(jié)構(gòu)不同迎角下橫向變形曲線

對比圖9、圖10數(shù)據(jù)可知，當(dāng)舵翼結(jié)構(gòu)迎角較小時，以縱向變形為主，橫向變形較小，隨著迎角增大，橫向變形明顯增大并超過縱向變形。為了進(jìn)一步分析上述變形規(guī)律產(chǎn)生的原因，對氣動熱、氣動力以及耦合作用下的橫向、縱向變形進(jìn)行分析。

舵翼結(jié)構(gòu)縱向變形如圖11所示，單純氣動熱作用下的變形量略大于耦合作用，二者變化規(guī)律一致，相差很?。欢鴼鈩恿Χ嬉斫Y(jié)構(gòu)造成壓縮變形，變形量隨迎角增大而降低，由于變形量非常小，對整體縱向變形的影響可以忽略。因此，舵翼結(jié)構(gòu)縱向變形主要由氣動熱的作用引起，氣動力作用極小，氣動熱與氣動力的耦合作用不明顯。

圖11 縱向變形曲線

舵翼結(jié)構(gòu)橫向變形如圖12所示，在氣動熱因素作用下變形為0.400 42～0.459 53 mm；在氣動力作用下橫向變形量較縱向變形量明顯增大，但仍然很小，僅為0.000 01～0.000 08 mm。當(dāng)氣動熱與氣動力耦合作用時，變形量明顯增大，達(dá)到為0.484 30～0.633 66 mm，可見耦合效應(yīng)在變形量上體現(xiàn)得比較明顯。在單純氣動力作用下，舵翼構(gòu)件并沒有發(fā)生明顯的彎曲變形。在氣動熱耦合作用下，舵翼結(jié)構(gòu)溫度急劇升高，材料物性參數(shù)發(fā)生明顯變化，彈性模量大幅降低，導(dǎo)致發(fā)生了更大的變形。彎曲變形改變了結(jié)構(gòu)的氣動外形，使溫度場進(jìn)一步升高，一方面加劇了因舵翼結(jié)構(gòu)熱膨脹和熱彎曲導(dǎo)致的變形，另一方面繼續(xù)降低了材料的彈性模量，使結(jié)構(gòu)產(chǎn)生了更大的變形。在二者的耦合作用，橫向變形量較單獨氣動熱和單獨氣動力作用下有著明顯的升高，迎角越大，氣動力與氣動熱作用的耦合效果越明顯。

圖12 橫向變形曲線

3.3 溫度分析

對舵翼結(jié)構(gòu)施加耦合工況，圖13為不同迎角狀態(tài)下的溫度分布云圖，舵翼結(jié)構(gòu)前緣處氣動熱效果最明顯，溫度最高，沿尾部方向逐漸降低。由于舵翼結(jié)構(gòu)前緣尖角相對氣流方向不是對稱的，上下面迎風(fēng)角度不同，導(dǎo)致溫度分布也不相同。下表面與來流形成的角度更大，氣動熱更劇烈，溫度更高。隨著迎角增大，下表面氣動熱效果進(jìn)一步加劇，溫度進(jìn)一步升高。

圖13 舵翼結(jié)構(gòu)溫度云圖

對舵翼結(jié)構(gòu)溫度場受氣動熱力耦合效應(yīng)的影響規(guī)律進(jìn)行研究，不同耦合因素下最高溫度隨迎角變化曲線如圖14所示。在3種耦合狀態(tài)下，溫度都隨迎角增大而升高。在只考慮氣動熱因素時，整體溫度最低；在無結(jié)構(gòu)變形耦合因素中，只考慮了舵翼結(jié)構(gòu)的熱膨脹而不考慮氣動力導(dǎo)致的變形，這時溫度較只考慮氣動熱時略高，可見單純熱膨脹變形能夠小幅度增加氣動熱效應(yīng)；在完全耦合因素中，溫度明顯大于上述2個工況，且迎角越大差值越高。由于在大迎角姿態(tài)飛行時，舵翼結(jié)構(gòu)受氣動熱和氣動力作用發(fā)生變形，使迎角進(jìn)一步增大，導(dǎo)致氣動加熱作用更加劇烈，進(jìn)而升高了舵翼結(jié)構(gòu)的溫度。可見，迎角越大，氣動熱與氣動力耦合效應(yīng)對結(jié)構(gòu)溫度場的影響就越明顯。

圖14 舵翼結(jié)構(gòu)溫度曲線

4 結(jié)論

1) 在近空間飛行狀態(tài)下，舵翼結(jié)構(gòu)頭部尖角處的溫度、變形最大，隨著飛行迎角增大，溫度升高，縱向變形減小，橫向變形增大。

2) 在近空間飛行狀態(tài)下，舵翼結(jié)構(gòu)會發(fā)生縱向拉伸和橫向彎曲變形，縱向拉伸變形主要由氣動加熱引起，橫向變形主要由氣動力作用引起。

3) 在近空間飛行狀態(tài)下，舵翼結(jié)構(gòu)氣動熱、氣動力的耦合效應(yīng)使變形與溫度增大；飛行迎角越大，耦合效應(yīng)越明顯。因此，當(dāng)舵翼結(jié)構(gòu)以大迎角狀態(tài)飛行時，氣動熱力環(huán)境會更加惡劣，在熱防護(hù)設(shè)計和結(jié)構(gòu)設(shè)計時應(yīng)充分考慮耦合效應(yīng)。