基于特征正交分解的一類瞬態非線性熱傳導問題的新型快速分析方法1)

朱強華 楊 愷 梁 鈺 高效偉

?(南昌航空大學飛行器工程學院,南昌 330063)

?(大連理工大學航空航天學院,工業裝備結構分析國家重點實驗室學,大連 116024)

引言

瞬態非線性熱傳導問題在實際工程中普遍存在,這是因為材料的熱物性參數通常是隨溫度變化的.由于這類問題的復雜性,很難得到其解析解,一般采用有限差分法[1-2]、有限元法[3-5]、有限體積法[6]等傳統的數值方法得到其近似解.鑒于這類問題的重要性,研究者仍不斷致力于發展新的求解手段.比如,顧元憲等[7]采用精細積分法求解瞬態非線性熱傳導問題,推導了熱傳導方程精細積分的具體列式;Tang 等[8]采用半邊界法(half-boundary method,HBM)求解核燃料棒的一維瞬態非線性熱傳導問題;Yang等[9]采用徑向積分邊界元法(radial integration boundary element method,RIBEM)求解導熱系數隨溫度變化形成的瞬態非線性熱傳導問題;Cui 等[10]采用Gao 等[11]近年提出的單元微分法(element differential method,EDM)求解多維瞬態非線性熱傳導問題.然而,這些數值方法在對工程中的大型、復雜系統進行求解時,得到的離散方程規模巨大,其自由度數量可達數百萬甚至千萬,這一方面對計算機的性能和存儲容量提出了極高的要求,另一方面需要耗費大量的計算時間,無法滿足一些需要對溫度進行實時控制或快速預測的領域的要求.為此,特征正交分解(proper orthogonal decomposition,POD)作為一種有效的模型降階方法,能夠將高維系統轉變為低維系統,從而大大縮短計算所需要的時間,因此是解決這個問題的重要途徑.

POD 方法又被稱為主成分分析(principal component analysis,PCA)或奇異值分解(singular value decomposition,SVD),在流體力學[12-13]、結構動力學[14-16]、最優控制[17-18]等諸多領域都得到廣泛應用.在傳熱學領域,Biaecki 等[19]采用POD 方法對瞬態線性熱傳導問題的有限元模型進行降階,使計算效率得到顯著提高;Zhang 和Xiang[20]將POD 方法和無網格法結合建立了瞬態線性熱傳導問題的降階模型,在保持計算精度的同時使計算時間大為減少;Gao 等[21]提出了一種多域POD 方法,能夠對由多種常物性材料組成的復雜求解域的瞬態熱傳導過程進行快速分析;胡金秀和高效偉[22]建立了變系數瞬態熱傳導問題邊界元格式的POD 降階模型,解決了邊界元法計算速度慢、算法穩定性差的問題;馮俞楷等[23]提出從溫度場的梯度提取POD 模態建立瞬態線性熱傳導問題的降階模型,不僅提高了計算效率,還具有更強的外推能力.這些研究展示出基于POD的模型降階方法在處理線性熱傳導問題時的高效性及其與各種數值方法結合的廣泛適用性.

雖然POD 是依賴信號來重構系統的一種線性處理方法,但研究表明它對于非線性系統也有一定的應用價值[24-25].目前國內外對瞬態非線性熱傳導問題的模型降階研究尚非常缺乏.Fic 等[26]初步研究了對瞬態非線性熱傳導問題采用POD 方法建立其降階模型的可行性,發現無需重新計算POD 模態;Han 等[27]針對瞬態變物性熱傳導問題發展了基于適體坐標(body-fitted coordinate,BFC)和有限體積法的POD 降階模型.盡管降階分析能夠準確地預測非線性熱傳導過程,但是對計算效率的提升效果卻遠不及其對線性熱傳導問題的那樣顯著,通常要小1～2個數量級.其原因是在瞬態非線性熱傳導問題的迭代計算過程中,每個迭代步都必須更新降階模型的系數矩陣,其關聯計算會耗費大量時間.于是POD降階模型的計算效率成為制約其在瞬態非線性熱傳導問題中應用的關鍵,也成為了近年來的研究熱點.Gaonkar 和Kulkarni[28]采用兩級離散和多級格式相結合的方法(two level discretization-multilevel scheme,TLD-MS)有效提高了瞬態非線性熱傳導問題的POD降階模型計算效率,但該方法局限于結構化網格且存在網格匹配性要求以及插值引起額外的精度損失;Feng 等[29]采用組合近似方法(combined approximations,CA),而Ding 等[30]則采用等幾何分析和獨立系數相結合的方法(isogeometric analysis-independent coefficients,IGA-IC),分別建立了瞬態非線性熱傳導問題的非POD 降階模型,雖然提高了計算效率,但在全階模型的自由度數量較少時效果亦不明顯.

在前期的研究中[31],我們發現可以采用瞬態線性熱傳導分析中得到的POD 模態對相同求解域的瞬態非線性熱傳導問題建立降階模型進行近似計算,這雖然簡化并加快了降階模型的構建,但并未涉及到如何解決非線性降階模型計算耗時的關鍵問題.本文對此展開了深入系統的高效算法研究,結構安排如下:第二部分給出瞬態非線性熱傳導問題的有限元全階模型及其計算方法;第三部分先采用POD 方法建立瞬態非線性熱傳導問題的降階模型,并簡要介紹其常規計算方法,然后重點闡述新型計算方法的基本理論;第四部分通過二維和三維算例驗證本文提出的瞬態非線性熱傳導降階模型新型計算方法的準確性和高效性;最后進行總結.

1 瞬態非線性熱傳導方程的有限元法

1.1 瞬態非線性熱傳導問題的控制方程

考慮材料的導熱系數隨溫度變化形成的非線性熱傳導問題,在求解域? 內其非穩態、無內熱源時的控制方程可表示為

式中,ρ 為材料的密度(kg/m3),c為材料的比質量熱容(J/(kg·K)),T(x,t)表示t時刻空間坐標x處的溫度(K),k(T)表示材料隨溫度T變化的導熱系數(W/(m·K)),t0為初始時刻(s).

方程(1)是一類拋物型偏微分方程,其定解條件包括初值條件

和邊值條件

式中,T0為初始時刻的溫度場,Γ=Γ1+Γ2+Γ3為求解域? 的邊界,其中Γ1是第一類邊界(Dirichlet 條件),是在邊界Γ1上給定的溫度,Γ2是第二類邊界(Neumann 條件),是在邊界Γ2上給定的熱流密度(W/m2),Γ3是第三類邊界(Robin 條件),是在邊界Γ3上給定的對流換熱系數(W/(m2·K)),Ta為環境溫度,ni(i=1,2,3)是邊界Γ 外法線的方向余弦.當和不隨時間變化時稱為定常邊界條件,否則稱為時變邊界條件.

這里需要說明的是,由于本文提出的模型降階快速分析方法目前僅適用于解決第二類和第三類邊界條件的瞬態非線性熱傳導問題,因此本文暫時不涉及第一類邊界條件的相關問題.

1.2 瞬態非線性熱傳導方程的有限元離散格式

應用有限元法求解上述瞬態非線性熱傳導問題,方程(1)及其定解條件(2)和(3)的弱解形式為

式中,N為權函數.由于所構造的權函數需要滿足正交的特性,這在整個求解域內是難以實現的.為此有限元法通過將求解域離散成若干個微元子域(即網格單元),在各子域內構造相應的權函數以求得其溫度場的近似解

式中,Te(t)為單元e的近似溫度,n為單元e所含的節點數量,和分別為單元e內節點i的權函數和溫度,Ne和Te(t)分別為單元e內各節點的權函數向量和溫度向量,nelem為求解域? 的離散單元數量.通過Galerkin 投影可以得到單元e的非線性代數方程組

式中,Ce,Ke和Fe分別是單元e的熱容矩陣、熱傳導矩陣和溫度載荷向量,其計算式如下

顯然,Ce和Ke都是n×n方陣,而Fe是n×1 列陣.

于是,上述瞬態非線性熱傳導問題的有限元離散格式可由形如式(6)的求解域? 內各單元的非線性代數方程組組集而成,即

式中,T(t)為t時刻的節點溫度向量

其中,nnode為求解域? 內的節點總量,而C,K和F分別是熱傳導系統的總體熱容矩陣、總體熱傳導矩陣和總體溫度載荷向量

其中,C和K都是nnode×nnode方陣,而F是nnode×1列陣.

由式(8)可以解得求解域? 的溫度場.由于瞬態非線性熱傳導控制方程(1)的有限元模型的方程數量等于節點總量(即自由度數量),當求解域? 的規模較大、形狀復雜時其內部須劃分數量龐大的網格單元,導致求解需要很大的存儲資源并耗費大量的計算時間.

1.3 瞬態非線性熱傳導有限元模型的計算方法

式(8)中總體熱傳導矩陣K是與溫度T(t)相關的,因此該式是一個非線性微分方程組.本文采用無條件穩定的Euler 隱式格式

將其轉化為非線性代數方程組,并應用N-R 方法進行迭代求解.

已知時刻tn的溫度場Tn,第i次迭代后的殘值可表示成

將殘值R應用泰勒展開,令有

式中,KT為切向剛度矩陣,KT=?R/(?T),于是

當修正溫度?T或殘值R達到收斂標準時,時刻tn+1的溫度場為

2 瞬態非線性熱傳導方程的模型降階法

2.1 瞬態非線性熱傳導問題的降階模型

對于需要實時控制或快速計算的工程領域,有限元等常規的數值方法不能滿足要求.因為實際問題的自由度數量高達上千萬甚至更多,求解這種超大規模的非線性方程組需要很長的時間.為此需要建立降階模型以減小計算量,從而達到實時控制或快速計算的目標.根據本課題組前期的研究[31],基于特征正交分解方法建立的降階模型能夠對瞬態非線性熱傳導過程進行準確預測,但是由于非線性問題的特殊性,降階模型在計算效率的提升方面不夠顯著,本文致力于解決這一問題.為保持文章的完整性,下面簡要介紹針對瞬態非線性熱傳導問題建立的降階模型,而關于特征正交分解方法的詳細內容可參閱文獻[32].

首先,通過實驗測量、數值模擬等方法獲得瞬態非線性熱傳導問題在若干個時刻點上的溫度場,通常稱其為瞬像,將所有瞬像按照時間順序排列成瞬像矩陣A=[T(t1)T(t2)···T(tl)]nnode×l.

其次,對瞬像矩陣A進行特征正交分解,獲得一組正交基矢量,通常稱其為POD 模態,截取前r階能夠捕捉絕大部分能量的模態構成模態矩陣Φ=

最后,利用模態矩陣Φ可以將任意時刻的溫度場T(t)表示成

2.2 瞬態非線性熱傳導問題降階模型的常規計算方法

采用隱式時間推進方法建立的瞬態非線性熱傳導降階模型的離散格式可表示成

然而需要強調的是,由式(16)可見此降階模型是不完備的,即仍然取決于原物理空間中的高維溫度向量Tn+1,而不能直接通過當前時刻得到的低維向量來計算,使得在式(18)的計算過程當中每個迭代步都要出現下列額外的計算環節:

(2)根據新的溫度場更新導熱系數,利用式(7b)和式(10b)重新計算各單元的熱傳導矩陣并組集成非線性熱傳導系統的總體熱傳導矩陣Kn+1.

(3)通過式(17b)的矩陣運算將Kn+1轉換成非線性熱傳導降階模型的系數矩陣

由于有限元全階模型的自由度數量較大,這些額外的數據轉換、矩陣組集和矩陣運算將耗費不少的計算時間,使得降低模型的小規模非線性代數方程組的求解帶來的時間收益被抵消.在后面的算例驗證部分我們可以看到,在有限元全階模型的規模不太大時,采用常規方法求解降階模型對計算效率的提升效果并不理想.

2.3 瞬態非線性熱傳導問題降階模型的新型計算方法

本小節將介紹瞬態非線性熱傳導降階模型的一種新型高效計算方法,它是直接通過對上述常規計算方法中兩個耗時環節的技術處理以達到提升計算效率的目的:一方面,采用單元預轉換方法(element preconversation method,EPM)壓縮降階模型系數矩陣的計算時間;另一方面,采用多級線性化方法(multilevel linearization method,MLM)消除非線性方程組的迭代計算過程.這種新方法解決了導熱系數隨溫度變化情況下EPM 和MLM 的結合問題,以及如何對換熱邊界單元的熱傳導矩陣進行相應計算.下面對此作詳細闡述.

2.3.1 EPM 理論與計算

這要求在有限元全階模型的單元層面就預先將單元熱傳導矩陣Ke進行轉換,其計算式為

式中,Φe是模態矩陣Φ中單元e的對應部分,Φe是n×r矩陣.顯然,與一樣,也是r×r方陣.

將式(7b)改寫成

于是,式(22)可以簡化成

將上式依次代入式(21)、式(20),可得

當瞬態非線性熱傳導問題的有限元全階模型和對流換熱邊界確定后,都可以預先計算并存儲備用.在迭代過程中,采用式(27)計算降階模型的系數矩陣時,只需要根據新的溫度場計算出各個單元的平均導熱系數即可.因此,與2.2 節的常規計算方法相比,EPM 能夠有效壓縮的計算時間,進而提高降階模型的計算效率.

2.3.2 MLM 理論與計算

MLM 的主要思想是通過插值技術取代非線性熱傳導系統的非線性項,得到其近似的線性方程組,從而利用前面若干級時刻點上已知的溫度場直接計算出當前時刻的溫度場.方便起見,將式(8)表示的瞬態非線性熱傳導問題的有限元全階模型改寫成

式中,系數陣W=CK(T),載荷向量f=CF(t).

本文采用無條件穩定的Dupont II 多級線性化格式來求解式(28),此格式為[33]

式中,W?表示插值溫度T?下的系數矩陣W,f?表示插值時間t?時的載荷向量f,T?和t?按下式計算

將式(15)代入式(29),然后方程兩邊同乘ΦT就可以得到瞬態非線性熱傳導降階模型的多級線性化格式

式中,Ir為r階單位矩陣,而分別為

由此降階模型的多級線性化格式可見,在前兩個時刻tn和tn+1的已知時,可先利用式(15)計算原物理空間的溫度場Tn和Tn+1,然后根據式(30)計算插值溫度T?和插值時間t?,繼而求出W?和f?,并利用式(32)將它們轉換為和就可以通過式(31)直接計算時刻tn+2的得到其溫度場Tn+2.

實際計算時,初始時刻t0的溫度場是已知的,而第二個時刻t1的溫度場未知,需要通過2.2 節的常規方法或者基于前面EPM 方法改進的常規方法求出,此后各個時刻的溫度場均可由多級線性化格式直接計算得到,消除了耗時的迭代過程,因而有利于提高降階模型的計算效率.

2.3.3 EPM 和MLM 的結合技術

鑒于EPM 和MLM 分別是針對常規方法中兩個不同環節做出的改進,這就為兩者結合以獲得降階模型更高的計算效率提供了可能.然而,其結合的困難在于對MLM 方法中式(32)的處理,將W?和f?代入其中可見

該式要求進行總體熱容矩陣的求逆和總體矩陣及其與模態矩陣之間的多重乘法運算,這與EPM 方法從單元層面的預處理方式是不協調的.

本文針對材料的比熱容為常數、導熱系數隨溫度變化這一類型的瞬態非線性熱傳導問題,提出有限元全階模型中的各個單元對降價模型系數矩陣亦有其相應的貢獻,貢獻大小用表示,即

由式(7a)得到的單元熱容矩陣Ce為協調質量熱容矩陣,為了計算式(35)中的單元矩陣須將其轉化為集中質量熱容矩陣,這樣由式(10a)得到的總體熱容矩陣C為對角矩陣,可寫為

由于比熱容為常數,于是C是不隨溫度變化的恒常矩陣,其逆矩陣亦為對角矩陣,此時可由下式進行計算

將式(24)代入式(37),并補充單元處于對流換熱邊界上時存在的單元熱傳導矩陣修正項可得

式中,重復指標p表示執行愛因斯坦求和約定,對二維問題和三維問題p分別取2 和3.

將上式依次代入式(35)、式(34),可得

這樣,就將EPM 和MLM 結合起來了.在用式(31)降階模型多級線性化格式進行計算時,只需要根據前兩個時刻已經得到的溫度場計算出插值溫度和插值時間即可,其余步驟和過程是與EPM 完全相同的.由此可見,這種針對瞬態非線性熱傳導降階模型發展的新型EPM-MLM 復合方法,與EPM 具有良好的一致性,其區別僅在于將和替換成了和簡單易行,便于編制規范統一的程序和軟件.

3 算例分析

本部分通過二維和三維瞬態非線性熱傳導算例的計算和分析,驗證本文所提出的降階模型新型計算方法的可行性、準確性和高效性.為方便對比說明,下面將瞬態非線性熱傳導降階模型的常規計算方法和新型計算方法分別命名為POD1 和POD2.

3.1 算例1:方板的瞬態非線性熱傳導問題

考慮如圖1 所示邊長為0.5 m 的二維方板,板材的密度ρ=7850 kg/m3、比熱容c=460 J/(kg·?C),其導熱系數隨溫度線性變化,k=65-0.05TW/(m·?C).方板的初始溫度為20?C,其右側為對流換熱邊界,換熱系數125 W/(m2·?C)、環境溫度1120?C;左側為放熱邊界,熱流密度10 000 W/m2;上、下側均為絕熱邊界.

圖1 方板的幾何模型及熱邊界條件Fig.1 Geometric model and thermal boundary conditions of the square plate

方板的有限元模型如圖2 所示.計算域采用線性三角形單元進行離散,各側邊均布21 個節點,總計888 個單元、485 個節點.為便于后處理分析,取方板的左側中點A、中心點B 和右側中點C 等三個點作為參考點.

為建立方板瞬態非線性熱傳導問題的降階模型,本文采取數值模擬的方法構造瞬像矩陣:時間步長取為50 s,通過有限元法計算0～20 000 s 時間范圍內方板在上述定常熱邊界條件下溫度場的變化,選取t=1000,2000,···,20 000 s 共計20 個時刻點上方板的溫度場構成瞬像矩陣A.對A進行特征正交分解分析,截取前4 個特征值(即r=4)對應的POD 模態構造模態矩陣Φ.利用式(15)可以得到僅包含4 個非線性方程的降階模型,相比于原來由485 個非線性方程組成的有限元全階模型,方程組的規模顯著減小.

圖2 方板的有限元模型Fig.2 Finite element model of the square plate

圖3 給出了有限元全階分析、常規降階分析和新型降階分析三種方法得到的方板在10 000 s 時刻的溫度場對比,可以看到三者是高度一致的.圖4 給出了方板上3 個參考點的溫度時變曲線,它將常規降階分析和新型降階分析兩種方法的計算結果與商業有限元軟件ANSYS 對此全階模型的計算結果(圖中以實線表示)進行比較,可見新型降階算法POD2 與常規的降階算法POD1 一樣,在所研究的整個時間范圍內能夠得到與有限元方法吻合的結果.因此,這些研究結果表明基于POD 的模型降階方法能夠用于瞬態非線性熱傳導問題的分析并得出與有限元方法相同的結果,而所提出的新型算法POD2 也可以得到正確的結果.

圖3 三種算法10 000 s 時刻方板的溫度場對比Fig.3 Comparison of temperature distributions in the square plate at t=10 000 s using different algorithms

圖3 三種算法10 000 s 時刻方板的溫度場對比(續)Fig.3 Comparison of temperature distributions in the square plate at t=10 000 s using different algorithms(continued)

圖4 參考點處溫度隨時間的變化曲線Fig.4 Temperature evolutions at reference points

下面進一步對降階模型新型算法的計算精度進行定量分析.首先,在方板的上邊線上取等間距的6個點,將其10 000 s 時刻的計算溫度和相對誤差列于表1.相對誤差是以有限元全階模型的計算結果為基準的誤差,其計算式為

可見,該時刻下POD2 算法的相對誤差相比于POD1算法略微增加,這是因為POD2 算法不僅包含POD1算法的模態截斷誤差,還附加了線性化誤差.但是POD2 算法的相對誤差不超過0.2%,仍然保持了比較高的精度,因而這種新型算法的計算結果是準確的.

表1 方板上邊線的計算溫度及相對誤差Table 1 Temperatures and relative errors on the top side of the square plate

其次,考慮所研究時間范圍內各個時刻的方板整體溫度的降階模型計算誤差,采用均方根誤差(root mean square error,RMSE)這一指標來衡量其大小

圖5 給出了RMSE 的時變曲線,可見降階模型的計算誤差在經過起始很短時間的脈動后迅速趨近于零,其峰值僅為1.1%.圖中也顯示出新型算法POD2 的誤差略大于常規算法POD1,其原因如前所述.

圖6 對有限元全階分析、常規降階分析和新型降階分析3 種方法花費的計算時間進行了比較,可以直觀看到POD1 用的計算時間較FEM 略有減少,而POD2 用的計算時間則大幅減少,這表明采用常規算法分析瞬態非線性熱傳導降階模型時的計算效率并不比有限元全階分析有明顯改進,而采用新型算法對降階分析計算效率的提高是非常顯著的.

圖5 方板溫度的均方根誤差隨時間的變化曲線Fig.5 Temporal variation of the RMSE for Example 1

圖6 各種算法所用計算時間的比較Fig.6 Comparison of computational times cost by various algorithms for Example 1

下面改變方板計算域的節點密度,使有限元全階模型的自由度數量從129 逐漸增加到2941,分析這個因素對降階模型POD1 和POD2 兩種算法計算效率的影響.表2 中列出了不同自由度下各種算法所花費的計算時間,同時給出了POD1 和POD2 的倍率數值.倍率是有限元全階分析與POD 降階分析所用計算時間的比值,其大小可以衡量降階分析對計算效率提升效果的強弱.由表可見,常規算法POD1 的計算效率提高得不多,尤其是在自由度數量較小時,降階分析所用的計算時間幾乎與全階分析的相近.與之相比,新型算法POD2 的計算效率則提高了2 個數量級,而且在自由度數量較小時,降階分析所用的計算時間也是大幅減少.此外,隨著自由度數量的增大,新型算法POD2 對計算效率的增幅愈加顯著.

表2 不同自由度時各種算法所用的計算時間Table 2 Computational times cost by various algorithms under different DOFs

最后,考慮到非線性熱傳導問題的有限元全階分析往往要花費較長的計算時間,而且實際工程上面對著各種復雜的熱邊界條件,測試也比較的困難,導致獲取相應問題的模態矩陣建立其降階模型既耗時又不易.如果能夠用簡單的定常熱邊界條件下建立的降階模型對相同計算域在復雜時變熱邊界條件下的瞬態非線性熱傳導過程進行分析和預測,將是非常有意義和應用價值的.為此,下面采用圖1 中定常熱邊界條件下建立的方板瞬態非線性熱傳導降階模型來分析其在新的時變熱邊界條件下的熱傳導過程,以檢驗這種方法的可行性.

保持方板右側對流換熱條件和上、下側絕熱條件不變,將其左側改為時變的放熱熱流邊界,如圖7所示,有二種情況:

(a)時變邊界1:熱流密度按正弦波

的形式光滑連續變化,在0～20 000 s 的時間范圍內周期性變化4 次.

(b)時變邊界2:熱流密度按三角波的形式非光滑連續變化,波動幅度增大為時變邊界1 的兩倍,而在0～20 000 s 的時間范圍內周期性變化次數則減少到兩次.

采用新型算法POD2 對前面定常熱邊界條件下建立的方板瞬態非線性熱傳導降階模型進行計算,得到兩種新的時變邊界條件下方板的溫度場,圖8 所示為3 個參考點上的溫度變化.由圖可見,點A 處的溫度明顯呈現出波動狀升溫的過程,因為該點就處于左側時變熱邊界上,甚至在時變邊界2 的條件下,處于方板右側邊界上的點C 也由于左側大幅度變化的熱流邊界而使其溫度呈現出輕微的波動.為驗證該降階模型計算結果的準確性,圖中同時以曲線形式給出了ANSYS 對相應邊界條件下方板溫度場的模擬結果,可見兩者是完全吻合的,因此用定常熱邊界條件下建立的瞬態非線性熱傳導降階模型能夠準確有效的對相同計算域在各種復雜的時變熱邊界條件下的熱傳導過程進行分析與預測.

圖7 方板左側施加的時變熱流Fig.7 Time-varying heat flux on the left side of the square plate

圖8 時變邊界下參考點處溫度的時變曲線Fig.8 Temperature evolutions at reference points under time-varying boundary conditions

3.2 算例2:散熱器的瞬態非線性熱傳導問題

考慮一個三維結構的散熱器,如圖9 所示,其底面為邊長0.2 m 的方形,高為0.15 m.3 個翅片的厚度及其間距均為δ,δ=0.04 m,翅片高為0.1 m.材料的密度ρ=7800 kg/m3、比熱容c=400 J/(kg·?C)、導熱系數隨溫度線性變化,k=60+0.05TW/(m·?C).散熱器的初始溫度為20?C,其底面是與發熱體接觸的加熱面,熱流密度50 000 W/m2;翅片與環境氣體間進行對流換熱,換熱系數100 W/(m2·?C)、環境溫度20?C;其余表面均為絕熱面.

圖9 散熱器的幾何模型及熱邊界條件Fig.9 Geometric model and thermal boundary conditions of the radiator

散熱器的有限元模型如圖10 所示.計算域采用六面體單元進行離散,總計4400 個單元、5796 個節點.為便于后處理,取中部翅片的角點A和B作為參考點,以及邊線CD進行溫度分析.

圖10 散熱器的有限元模型Fig.10 Finite element model of the radiator

與算例1 一樣,為建立散熱器的瞬態非線性熱傳導降階模型,采取數值模擬的方法構造瞬像矩陣:時間步長取5 s,通過有限元法計算0～500 s 時間范圍內散熱器溫度場的變化,選取t=25,50,···,500 s 共計20 個時刻點上的溫度場構成瞬像矩陣A.對A進行特征正交分解分析,截取前5 個特征值(即r=5)對應的POD 模態構造模態矩陣Φ.利用式(15)可以得到僅包含5 個非線性方程的降階模型,相比于由5796 個非線性方程組成的有限元全階模型,方程組的規模急劇減小.

圖11 給出了有限元全階分析、常規降階分析和新型降階分析3 種方法得到的散熱器在末時刻500 s 時的溫度場對比,可以看到三者基本上是相同的.圖12 給出了散熱器上A,B二個參考點的溫度時變曲線,圖中將POD1 和POD2 兩種降階模型算法的計算結果與商業有限元軟件ANSYS 對全階模型的計算結果進行了比較,同樣可以看到新型算法POD2與常規算法POD1 一樣能夠得到與有限元方法吻合的結果.因此本文提出的新型算法POD2 可以對瞬態非線性熱傳導問題進行準確的分析.

圖11 三種算法500 s 時刻散熱器的溫度場對比Fig.11 Comparison of temperature distributions in the radiator at t=500 s using different algorithms

圖12 參考點處溫度隨時間的變化曲線Fig.12 Temperature evolutions at reference points

在散熱器的邊線CD 上取等間距的6 個點,將其500 s 時刻的計算溫度和相對誤差列于表3 中.由表可見,與算例1 類似,該時刻下POD2 算法的相對誤差相比于POD1 算法略微增大.常規算法POD1 的最大相對誤差約為1.12%,而新型算法POD2 的最大相對誤差約為1.81%,因而這種新型算法仍然保持了比較高的精度,其計算結果是準確的.

表3 散熱器邊線CD 上的計算溫度及相對誤差Table 3 Temperatures and relative errors on the line CD of the radiator

圖13 給出了用POD1 和POD2 兩種降階模型算法計算得到的散熱器整體溫度的均方根誤差RMSE的時變曲線,可以看到降階模型的計算誤差經過起始時間的脈動,在約200 s 后緩慢穩定下降.POD1和POD2 的脈動峰值分別為0.01%和0.03%.降階模型的計算誤差在穩定段基本小于0.01%,新型算法POD2 的誤差略大于常規算法POD1.

圖14 對散熱器算例采用有限元全階分析、常規降階分析和新型降階分析3 種方法花費的計算時間進行了比較.由圖可見,由于該算例的自由度數量比較大,有限元全階分析需花費長達約862 s 的計算時間,此時常規降階分析和新型降階分析花費的計算時間均有不同程度減少,尤其是后者,用時僅為7 s,計算效率提高了將近123 倍.

圖13 散熱器溫度的均方根誤差隨時間的變化曲線Fig.13 Temporal variation of the RMSE for Example 2

圖14 各種算法所用計算時間的比較Fig.14 Comparison of computational times cost by various algorithms for Example 2

3.3 新型算法的特性分析

通過前面兩個算例的計算分析,驗證了本文提出的降階模型新型算法是準確有效的,特別是在計算效率方面,較常規算法有大幅度的提升.

與近年來發展的針對瞬態非線性熱傳導問題的其它一些降階算法相比,本文的新型算法也具有更高的計算效率.表4 中列有IGA-IC[30]、TLD-MS[28]和CA[29]等一些方法對自由度數量與本文算例盡量接近的算例所用的計算時間,二維問題和三維問題分欄給出.須說明的是,CA 方法的計算時間文獻中未直接給出,是從其圖上提取的,故標為約數.考慮到各個研究者所用計算機性能的不同,對計算時間是有影響的,故表中同時給出了倍率值(全階與降階計算時間的比值),以更客觀的評價各種方法對降階模型計算效率的提升效果.由表可見,本文的POD2 方法在計算效率上較其它方法是占優的,特別是三維問題.而二維問題中TLD-MS 方法的計算效率雖然接近POD2 方法,但該方法需要建立疏、密兩套網格,算法繁雜,既存在網格匹配的要求,又會因不同網格間的插值造成精度損失,不利于通用軟件的編制和實際工程的應用.

表4 新型算法和文獻中其他方法的計算效率對比Table 4 Comparison of computational efficiency among POD2 and other methods in literatures

4 結論

本文針對導熱系數隨溫度變化的一類瞬態非線性熱傳導問題,發展并提出了一種基于POD 模型降階技術的新型高效快速分析方法.通過單元預轉換方法和多級線性化方法的有效結合解決了非線性降階模型的常規計算方法在系數矩陣更新和迭代等計算環節的費時問題,從而能夠獲得很高的計算效率.

通過二維和三維算例分析,演示了瞬態非線性熱傳導降階模型新型計算方法的應用細節,通過與全階模型有限元法計算結果的比較,驗證了該新型算法的準確性,并通過與降階模型常規計算方法以及近年來發展出的其他新方法在計算時間上的比較,驗證了該新型算法的高效性.

研究結果表明,該新型算法具有良好的特性:其一,計算穩定性好,算法基于的隱式時間推進格式和多級線性化格式均為無條件穩定;其二,計算速度快,常規算法在全階模型自由度數量較小時計算時間沒有明顯減少,而新型算法無論自由度數量是大是小都能夠加快計算速度,并且自由度數量越多計算效率的提升效果越好,在本文自由度數量不超過6000 的小規模計算中其計算效率能提高2～3 個數量級;其三,計算通用性好,便于編制規范、統一的程序和軟件,有利于實際應用.

此外,采用常數熱邊界條件下得到的POD 模態建立的瞬態非線性熱傳導降階模型可以對相同計算域在各種復雜的光滑/非光滑時變熱邊界條件下的非線性熱傳導過程進行快速、準確的分析和預測,這進一步拓展了本文所提降階模型新型算法在工程應用中的價值.