基于分區(qū)耦合的高超聲速熱防護結(jié)構(gòu)耦合傳熱分析

基于分區(qū)耦合的高超聲速熱防護結(jié)構(gòu)耦合傳熱分析

由于高速汽機和遠程機動性飛行在惡劣的加熱環(huán)境中的長期飛行，這將給設(shè)計師帶來許多新的問題和挑戰(zhàn)，如熱防護系統(tǒng)的內(nèi)部空間安排系統(tǒng)和熱氣動彈性的動態(tài)彈性。對于這些氣動加熱與結(jié)構(gòu)熱響應(yīng)強烈耦合的問題,傳統(tǒng)的熱環(huán)境和熱響應(yīng)分割獨立計算的方法已不再適用,必須采用流-固-熱多場耦合計算方法進行研究。多場耦合計算涉及各物理場的準確求解。目前以ANSYS和NASTRAN為代表的非線性有限元分析技術(shù)已較完善,可以分析包含各種幾何與材料非線性、輻射、相變和復(fù)雜結(jié)構(gòu)等結(jié)構(gòu)動力學(xué)及傳熱問題。對于流場求解,基于Navier-Stokes方程的計算流體力學(xué)(CFD)軟件也較為成熟,如FLUENT和CFD-FASTRAN等,可以模擬從低速到高速以及包含化學(xué)反應(yīng)的非平衡流等復(fù)雜流動問題。因此,開發(fā)通用的耦合計算平臺,充分利用現(xiàn)有軟件的強大功能,可以為流-固-熱多場耦合的復(fù)雜工程問題提供有效的解決手段。耦合計算平臺的關(guān)鍵問題包括以下兩個方面:①載荷傳遞。對于位移和溫度等場變量一般采用多項式插值、最近鄰居(NearestNeighbor)以及樣條插值等方法計算,其算法較為成熟。對于壓力和熱流通量等載荷則需要滿足守恒性,即需要采用守恒插值(ConservativeRemapping)方法。②各場軟件之間的數(shù)據(jù)傳遞和同步。實現(xiàn)多個軟件之間的交互,需要從操作系統(tǒng)入手進行內(nèi)核級別的程序開發(fā)。對于ANSYS和Fluent等提供了源代碼級別二次開發(fā)功能的軟件實現(xiàn)起來較為容易,對于CFD-FASTRAN等未提供二次開發(fā)功能的軟件,則需要對操作系統(tǒng)有深入的了解。在高超聲速熱防護結(jié)構(gòu)(TPS)的耦合傳熱計算方面,由于實際結(jié)構(gòu)及傳熱過程的復(fù)雜性,工程上多采用分割獨立求解方法來解決,也有采用工程算法或整場方法來進行處理的。需要指出的是,在高溫條件下,材料的發(fā)射率非線性以及結(jié)構(gòu)尺寸等因素對傳熱過程均將產(chǎn)生顯著影響,但目前尚未見到相關(guān)報導(dǎo)。針對多場耦合計算平臺的關(guān)鍵問題,本文提出一種滿足局部守恒的載荷插值方法,根據(jù)已有的有限元和CFD分析軟件的不同特性建立統(tǒng)一的多場耦合計算平臺,并進行算例驗證。在此基礎(chǔ)上,對高超聲速巡航狀態(tài)下的類X-34頭部熱防護結(jié)構(gòu)進行耦合傳熱計算,并進一步探討防護層結(jié)構(gòu)尺寸及材料非線性對于熱防護性能的影響。1算法的基本思想對于流固耦合傳熱問題,在流固交界面上需要滿足溫度和熱流的一致性,但二者在計算前對流體和固體都是未知的,因此需要采用迭代方法求解。目前普遍采用的方法是先假設(shè)適當?shù)谋诿鏈囟?將此壁面溫度作為流體邊界條件求解得到新的壁面熱流,將新的熱流代入固體傳熱方程求解得到新的邊界溫度,反復(fù)循環(huán)直到迭代達到收斂要求,從而得到最終解。根據(jù)耦合方式的不同可分為松耦合和緊耦合兩種。緊耦合需要在每個時間步上進行反復(fù)迭代計算來滿足邊界溫度及通量的一致性,計算量較大;相對而言,松耦合在一個時間步上僅計算一次,在時間步長足夠小的情況下其計算結(jié)果是可以接受的。流固耦合傳熱問題按照時間精度要求又可分為定常和非定常兩種,對于定常耦合問題需要采用緊耦合方式,非定常問題可選擇松耦合和緊耦合兩種。本文所述方法對二者均適用。耦合迭代的收斂準則通常采用慣用的均方根相對殘值公式,即√Ν∑i=1(f(n+1)i-f(n)i)2Ν∑i=1(f(1)i-f(0)i)2≤εtol(1)∑i=1N(f(n+1)i?f(n)i)2∑i=1N(f(1)i?f(0)i)2???????????ue001?ue000ue000ue000ue000≤εtol(1)式中:fi可取邊界溫度或熱流密度,下標i表示第i個流固耦合邊界節(jié)點;εtol為給定收斂誤差;N為邊界節(jié)點總數(shù);上標0,1,…,n+1表示流固耦合迭代步次數(shù)。由于熱流密度定義為溫度梯度的線性函數(shù),因此二者的收斂性是一致的?？紤]到一般結(jié)構(gòu)實際溫度范圍在102～103K量級,當溫度相對殘值達到10-5量級時溫度誤差的最大值在10-2～10-3K量級,因此收斂誤差取小于10-5的值即可。2單元特性fvm一般情況下,流體網(wǎng)格和結(jié)構(gòu)網(wǎng)格在邊界處不是一一對應(yīng)的,確保載荷的精確傳遞是保證耦合計算精度的關(guān)鍵。本文從線性結(jié)構(gòu)有限元法(FEM)和格心式有限體積法(FVM)的單元特性出發(fā),提出滿足局部守恒的通量插值方法,可以精確實現(xiàn)熱流通量的守恒傳遞。對溫度及位移按照流體網(wǎng)格邊界面的面心坐標位置,采用有限單元形函數(shù)進行插值,確保流固邊界溫度及位移的一致性,插值過程較為簡單。壓力和熱流屬于分布式載荷,其插值方法需要滿足守恒性要求,二者的計算方法是相同的,下面僅以熱流通量插值方法為例給出推導(dǎo)過程。2.1求解域的熱流通量的計算圖1為三維流固邊界上的一個結(jié)構(gòu)FEM邊界面和流體FVM邊界面的對應(yīng)關(guān)系示意圖,其中背景網(wǎng)格為流體網(wǎng)格。邊界上的熱流通量應(yīng)滿足如下守恒條件:(1)局部守恒性,即對單個單元滿足∫SkFkds=∫Skfkds(2)∫SkFkds=∫Skfkds(2)式中:Sk為有限單元k的外表面;Fk為Sk上的熱流通量;fk為有限體積網(wǎng)格對應(yīng)于Sk上的熱流通量。(2)整體守恒性,對整個求解域,滿足∫SFds=∫Sfds(3)∫SFds=∫Sfds(3)式中:S為整個流固耦合邊界;F為求解域結(jié)構(gòu)邊界的熱流通量;f為求解域流體邊界的熱流通量。在求解域完全重合的條件下,如果條件(1)能夠精確滿足,則條件(2)是自然滿足的。對于采用格心式FVM網(wǎng)格,式(2)可寫為∫SkFkds=Ν′∑j=1fk,jSk,j(4)式中:fk,j和Sk,j為與有限單元k相交的第j個有限體積網(wǎng)格的邊界熱流通量和面積;N′為與有限單元k相交的有限體積網(wǎng)格個數(shù)。2.2熱流通量的局部性質(zhì)以四邊形有限單元為例,單元溫度分布滿足Τ(s,t)=4∑i=1ΝiΤi(5)式中:Ti為節(jié)點i的溫度;s∈[-1,1]和t∈[-1,1]為自然坐標;Ni為形函數(shù),且Νi=14(1+sis)(1+tit)式中:si和ti為節(jié)點i的局部坐標。依據(jù)熱流通量的定義,熱流通量在單元上的分布滿足F=ks?Τ?s+kt?Τ?t(6)式中:ks和kt為局部坐標系中的導(dǎo)熱系數(shù)。假設(shè)t=-1邊為流固耦合邊界,顯然,熱流通量在邊界上滿足線性分布,不妨設(shè)為F(s)=as+b(7)將式(7)代入式(4)得到Sk∫1-1Fds=Ν′∑j=1fk,jSk,j積分后得到b=Ν′∑j=1fk,jSk,j2Sk又由fk滿足式(7),可得fk,1=as1+bfk,2=as2+b?fk,Ν′=asΝ′+b}(8)問題轉(zhuǎn)換為線性最小二乘問題,其求解方法較為簡單,不在此給出。由于三角形單元可以作為退化的四邊形單元對待,因此上述方法對于三角形單元同樣適用。2.3流固耦合邊界面心計算以8節(jié)點六面體單元為例,單元溫度分布滿足Τ(s,t,r)=8∑i=1ΝiΤi(9)式中:Νi=18(1+sis)(1+tit)(1+rir);s∈[-1,1]、t∈[-1,1]和r∈[-1,1]為自然坐標,si、ti和ri為節(jié)點i的局部坐標。單元熱流通量分布函數(shù)為F(s,t,r)=ks?Τ?s+kt?Τ?t+kr?Τ?r(10)假設(shè)r=-1面為流固耦合邊界,將式(9)代入式(10),可得該邊界面上熱流通量為二次函數(shù),記為F(s,t)=ast+bs+ct+d(11)將式(11)代入式(4)可得d=Ν′∑j=1fk,jSk,j4Sk與二維問題相同,最終轉(zhuǎn)換為線性最小二乘問題,即fk,1=as1t1+bs1+ct1+dfk,2=as2t2+bs2+ct2+d?fk,Ν′=asΝ′tΝ′+bsΝ′+ctΝ′+d}(12)式(12)中僅包含3個未知數(shù),可以通過3階矩陣的求逆公式精確計算。注意式(8)和式(12)都需要計算有限體積單元邊界面面心的局部坐標,對于非一對一的邊界網(wǎng)格需要查詢流體網(wǎng)格節(jié)點所對應(yīng)的固體網(wǎng)格單元?？紤]到邊界流體網(wǎng)格節(jié)點位移也是通過形函數(shù)插值得到的,因此相應(yīng)查詢計算在迭代過程中僅需進行一次,故可采用簡單的逐點查詢。本文采用Kdtree算法實現(xiàn)快速查詢。完成查詢后,將流體網(wǎng)格點投影到有限元單元計算相應(yīng)局部坐標,三維情況的重疊面積采用多邊形求交算法實現(xiàn)。3平臺構(gòu)建的關(guān)鍵技術(shù)傳統(tǒng)的基于多個軟件之間的耦合計算一般采用輸入文件方式,即:啟動軟件—讀取數(shù)據(jù)—計算—保存數(shù)據(jù)到磁盤。這種頻繁的軟件啟動和磁盤文件操作會極大降低計算速度。鑒于目前大多數(shù)結(jié)構(gòu)有限元分析和CFD軟件都提供了二次開發(fā)接口,如ANSYS的UPF和Fluent的UDF等,充分利用這些接口可以給耦合計算效率帶來質(zhì)的提升。本文通過自主開發(fā)的耦合計算平臺來協(xié)調(diào)各個軟件之間的數(shù)據(jù)傳遞和同步,實現(xiàn)多種軟件耦合計算。操作系統(tǒng)環(huán)境為WindowsXP,開發(fā)語言為C++。高效數(shù)據(jù)傳遞方式是解決耦合計算問題的關(guān)鍵之一,特別是對于非定常計算。平臺構(gòu)建的兩個關(guān)鍵技術(shù)如下所述:(1)數(shù)據(jù)傳遞。采用共享內(nèi)存方式實現(xiàn),對于提供了源代碼級別二次開發(fā)功能的軟件(如ANSYS和FLUENT等)可以直接在源程序中調(diào)用系統(tǒng)應(yīng)用程序編程接口(ApplicationProgrammingInterface,API)函數(shù)實現(xiàn)。對于不提供該功能的軟件(如FASTRAN),先通過系統(tǒng)鉤子(Hook)函數(shù)監(jiān)視磁盤文件操作,然后注入相同代碼實現(xiàn)數(shù)據(jù)傳遞。為了防止共享內(nèi)存訪問沖突,需要對數(shù)據(jù)訪問進行加鎖,采用互斥體(Mutex)實現(xiàn),典型的數(shù)據(jù)訪問流程如圖2所示。(2)同步方式。采用系統(tǒng)提供的互斥體和事件(Event)實現(xiàn)。由于平臺是構(gòu)建在各個軟件相同級別的層次上,因此即使相關(guān)軟件是并行運行的,共享內(nèi)存也是適用的。故本文方法適用于大規(guī)模問題并行計算。圖3所示為基于ANSYS和FLUENT的耦合傳熱計算流程圖。ANSYS和CFD-FASTRAN軟件的耦合流程與之類似,在此不做詳述。4計算4.1算例2:定常緊耦合+不同網(wǎng)格密度的湍流模型噴管幾何尺寸、流場以及壁面邊界條件參見文獻,噴管材料按照文獻取AISI302不銹鋼。耦合計算的結(jié)構(gòu)傳熱軟件為ANSYS,CFD軟件采用FLUENT,FEM單元選取4節(jié)點PLANE55單元?？紤]到噴管的軸對稱幾何特性以及外壁溫度較低(322K),噴管內(nèi)外壁均未考慮輻射作用。湍流模型選取SSTk-ω兩方程模型。此算例耦合計算方式為定常緊耦合。計算用流體網(wǎng)格為149×51結(jié)構(gòu)化網(wǎng)格。結(jié)構(gòu)有限元網(wǎng)格分別采用兩種不同網(wǎng)格密度以考核插值精度:①邊界流固節(jié)點一對一;②邊界結(jié)構(gòu)網(wǎng)格密度為流體網(wǎng)格的1/2。圖4所示為局部網(wǎng)格圖。圖5為兩種不同網(wǎng)格密度下的溫度云圖,兩種網(wǎng)格的溫度結(jié)果差別很小。耦合壁面上的溫度分布曲線及不同網(wǎng)格的壁面熱流通量相對殘值收斂曲線如圖6所示。顯然兩種網(wǎng)格的溫度分布和實驗結(jié)果非常接近。兩種網(wǎng)格的收斂速度差別不大,在20個迭代步即可收斂到10-7量級,隨后殘值收斂曲線的小幅震蕩是受限于FLUENT的收斂標準的緣故,此時FLUENT的解已收斂到指定量級10-6,如果降低FLUENT的收斂標準,迭代會繼續(xù)收斂下去。4.2實驗和元網(wǎng)格的計算本算例的結(jié)構(gòu)傳熱部分采用ANSYS軟件計算,高超聲速流場采用CFD-FASTARN軟件求解,迭代過程采用松耦合方法。結(jié)構(gòu)尺寸及來流參數(shù)參見文獻,來流壓力為648.1Pa,溫度為241.5K,馬赫數(shù)為6.47。結(jié)構(gòu)材料選取標準1Cr18Ni9Ti不銹鋼(即AISI321),考慮了材料的溫度非線性,有限元單元選取4節(jié)點PLANE13多場耦合單元。流場和結(jié)構(gòu)有限元網(wǎng)格如圖7所示,其中流體網(wǎng)格按照定常流場結(jié)果在激波位置進行了局部加密。湍流模型選取SSTk-ω兩方程模型,壁面第1層網(wǎng)格高度為1.0×10-5m。耦合時間步長為0.0001s,總計算時間為2s。非定常計算時,流場初始化為按照絕熱邊界計算的定常流場。結(jié)構(gòu)初始溫度為294.4K,分別在y=0m和x=0m邊界施加x和y方向位移約束。圖8為2s時刻的流體密度分布比較圖。圖中:ρ為流體密度;上半部分為實驗結(jié)果照片和文獻的計算結(jié)果。從圖中可以看出,本文激波位置與實驗結(jié)果較吻合。圖9為2s時刻的流場和結(jié)構(gòu)溫度分布云圖,其中結(jié)構(gòu)最大溫度出現(xiàn)在駐點位置,其值為444K,略低于實驗值465K。駐點熱流密度為6.656×105W/m2,與實驗值的6.7×105W/m2非常接近。圖10為壁面溫度和熱流分布曲線。圖中:T0為駐點溫度;θ為圓柱外表面與對稱面的夾角;q為壁面熱流密度;q0為駐點位置壁面熱流密度。從圖中可以看出,本文計算和實驗結(jié)果吻合較好。圖11為2s時刻的結(jié)構(gòu)熱應(yīng)力分布云圖,應(yīng)力準則采用vonMises理論。最大應(yīng)力出現(xiàn)在駐點位置,為442MPa,高于文獻的結(jié)果,這是由于本文的熱流密度結(jié)果更接近實驗值,高于文獻約38%,同時本文考慮了材料非線性作用的原因。4.3熱防護層厚度的影響對于長時間暴露在高超聲速流中的飛行器,其飛行時間已超過熱平衡所需時間,熱防護結(jié)構(gòu)需要采用內(nèi)部主動散熱系統(tǒng)。下面針對類X-34飛行器頭部金屬熱防護結(jié)構(gòu)進行定常耦合計算,考察其熱防護的性能。如圖12所示,熱防護結(jié)構(gòu)分為3層:最外層為76mm的LI-900防熱瓦,中間層為4.4mm應(yīng)變隔離墊(SIP),下層為1.6mm的鋁板(材料為Al2042)。材料的導(dǎo)熱系數(shù)、發(fā)射率均與壓力及溫度呈非線性關(guān)系,具體參數(shù)值參見文獻。本算例采用ANSYS和FASTARN進行耦合計算,考慮了材料非線性以及輻射效應(yīng)。結(jié)構(gòu)部分采用SHELL131層合板單元用于模擬熱傳導(dǎo),采用SUF152表面效應(yīng)單元模擬熱輻射。最內(nèi)層鋁板內(nèi)壁設(shè)置323K溫度約束。來流條件:35km高空大氣;迎角為23°;飛行馬赫數(shù)為6.3。結(jié)構(gòu)網(wǎng)格以及流體網(wǎng)格局部圖如圖13所示。流場壓力及壁面熱流密度分布如圖14所示,圖中:p為流場壓力。從圖中可以看出,最高壓力及熱流密度均位于頭部駐點區(qū)域,壓力值為2.925×104Pa,熱流密度值為1.62×105W/m2。底部迎風(fēng)面為高熱流區(qū)域,熱流密度在5×104W/m2量級上。圖15為壁面外表面以及LI-900防熱瓦底部的溫度分布云圖,駐點溫度為1119K。LI-900防熱瓦底部的最高溫度為489K,基體鋁板外表面溫度最高為423.1K,應(yīng)變隔離墊承受了50K的溫差,基體內(nèi)壁的熱流總值為4581W,即內(nèi)壁散熱系統(tǒng)需要4581W功率來確保結(jié)構(gòu)溫度在323K范圍內(nèi)。圖16為壁面熱流密度總殘值和流場密度相對殘值的迭代收斂曲線。ANSYS和CFD-FASTRAN的迭代次數(shù)為30次,CFD-FASTRAN總迭代步為1500,迭代頻率為50步。20個耦合迭代步時殘值降低到10-5量級。防護層厚度是熱防護結(jié)構(gòu)的主要參數(shù)之一,為考慮其影響,文中對20～120mm范圍內(nèi)的14種厚度情況進行了計算。如圖17為駐點溫度和內(nèi)壁冷卻系統(tǒng)功率隨防護層厚度的變化曲線。圖中:T0-NL表示非線性材料的駐點溫度;HF-NL表示非線性材料的駐點熱流。從圖中可以看出:駐點溫度隨厚度增加而增加,但總體變化在20K左右,在厚