朱強(qiáng)華 楊 愷 梁 鈺 高效偉

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

?(大連理工大學(xué)航空航天學(xué)院,工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室學(xué),大連 116024)

引言

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

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

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

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

1 瞬態(tài)非線性熱傳導(dǎo)方程的有限元法

1.1 瞬態(tài)非線性熱傳導(dǎo)問題的控制方程

考慮材料的導(dǎo)熱系數(shù)隨溫度變化形成的非線性熱傳導(dǎo)問題,在求解域? 內(nèi)其非穩(wěn)態(tài)、無內(nèi)熱源時(shí)的控制方程可表示為

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

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

和邊值條件

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

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

1.2 瞬態(tài)非線性熱傳導(dǎo)方程的有限元離散格式

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

式中,N為權(quán)函數(shù).由于所構(gòu)造的權(quán)函數(shù)需要滿足正交的特性,這在整個(gè)求解域內(nèi)是難以實(shí)現(xiàn)的.為此有限元法通過將求解域離散成若干個(gè)微元子域(即網(wǎng)格單元),在各子域內(nèi)構(gòu)造相應(yīng)的權(quán)函數(shù)以求得其溫度場的近似解

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

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

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

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

式中,T(t)為t時(shí)刻的節(jié)點(diǎn)溫度向量

其中,nnode為求解域? 內(nèi)的節(jié)點(diǎn)總量,而C,K和F分別是熱傳導(dǎo)系統(tǒng)的總體熱容矩陣、總體熱傳導(dǎo)矩陣和總體溫度載荷向量

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

由式(8)可以解得求解域? 的溫度場.由于瞬態(tài)非線性熱傳導(dǎo)控制方程(1)的有限元模型的方程數(shù)量等于節(jié)點(diǎn)總量(即自由度數(shù)量),當(dāng)求解域? 的規(guī)模較大、形狀復(fù)雜時(shí)其內(nèi)部須劃分?jǐn)?shù)量龐大的網(wǎng)格單元,導(dǎo)致求解需要很大的存儲資源并耗費(fèi)大量的計(jì)算時(shí)間.

1.3 瞬態(tài)非線性熱傳導(dǎo)有限元模型的計(jì)算方法

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

將其轉(zhuǎn)化為非線性代數(shù)方程組,并應(yīng)用N-R 方法進(jìn)行迭代求解.

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

將殘值R應(yīng)用泰勒展開,令有

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

當(dāng)修正溫度?T或殘值R達(dá)到收斂標(biāo)準(zhǔn)時(shí),時(shí)刻tn+1的溫度場為

2 瞬態(tài)非線性熱傳導(dǎo)方程的模型降階法

2.1 瞬態(tài)非線性熱傳導(dǎo)問題的降階模型

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

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

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

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

2.2 瞬態(tài)非線性熱傳導(dǎo)問題降階模型的常規(guī)計(jì)算方法

采用隱式時(shí)間推進(jìn)方法建立的瞬態(tài)非線性熱傳導(dǎo)降階模型的離散格式可表示成

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

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

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

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

2.3 瞬態(tài)非線性熱傳導(dǎo)問題降階模型的新型計(jì)算方法

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

2.3.1 EPM 理論與計(jì)算

這要求在有限元全階模型的單元層面就預(yù)先將單元熱傳導(dǎo)矩陣Ke進(jìn)行轉(zhuǎn)換,其計(jì)算式為

式中,Φe是模態(tài)矩陣Φ中單元e的對應(yīng)部分,Φe是n×r矩陣.顯然,與一樣,也是r×r方陣.

將式(7b)改寫成

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

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

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

2.3.2 MLM 理論與計(jì)算

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

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

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

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

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

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

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

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

2.3.3 EPM 和MLM 的結(jié)合技術(shù)

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

該式要求進(jìn)行總體熱容矩陣的求逆和總體矩陣及其與模態(tài)矩陣之間的多重乘法運(yùn)算,這與EPM 方法從單元層面的預(yù)處理方式是不協(xié)調(diào)的.

本文針對材料的比熱容為常數(shù)、導(dǎo)熱系數(shù)隨溫度變化這一類型的瞬態(tài)非線性熱傳導(dǎo)問題,提出有限元全階模型中的各個(gè)單元對降價(jià)模型系數(shù)矩陣亦有其相應(yīng)的貢獻(xiàn),貢獻(xiàn)大小用表示,即

由式(7a)得到的單元熱容矩陣Ce為協(xié)調(diào)質(zhì)量熱容矩陣,為了計(jì)算式(35)中的單元矩陣須將其轉(zhuǎn)化為集中質(zhì)量熱容矩陣,這樣由式(10a)得到的總體熱容矩陣C為對角矩陣,可寫為

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

將式(24)代入式(37),并補(bǔ)充單元處于對流換熱邊界上時(shí)存在的單元熱傳導(dǎo)矩陣修正項(xiàng)可得

式中,重復(fù)指標(biāo)p表示執(zhí)行愛因斯坦求和約定,對二維問題和三維問題p分別取2 和3.

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

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

3 算例分析

本部分通過二維和三維瞬態(tài)非線性熱傳導(dǎo)算例的計(jì)算和分析,驗(yàn)證本文所提出的降階模型新型計(jì)算方法的可行性、準(zhǔn)確性和高效性.為方便對比說明,下面將瞬態(tài)非線性熱傳導(dǎo)降階模型的常規(guī)計(jì)算方法和新型計(jì)算方法分別命名為POD1 和POD2.

3.1 算例1:方板的瞬態(tài)非線性熱傳導(dǎo)問題

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

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

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

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

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

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

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

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

圖4 參考點(diǎn)處溫度隨時(shí)間的變化曲線Fig.4 Temperature evolutions at reference points

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

可見,該時(shí)刻下POD2 算法的相對誤差相比于POD1算法略微增加,這是因?yàn)镻OD2 算法不僅包含POD1算法的模態(tài)截?cái)嗾`差,還附加了線性化誤差.但是POD2 算法的相對誤差不超過0.2%,仍然保持了比較高的精度,因而這種新型算法的計(jì)算結(jié)果是準(zhǔn)確的.

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

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

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

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

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

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

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

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

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

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

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

的形式光滑連續(xù)變化,在0～20 000 s 的時(shí)間范圍內(nèi)周期性變化4 次.

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

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

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

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

3.2 算例2:散熱器的瞬態(tài)非線性熱傳導(dǎo)問題

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

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

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

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

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

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

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

圖12 參考點(diǎn)處溫度隨時(shí)間的變化曲線Fig.12 Temperature evolutions at reference points

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

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

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

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

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

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

3.3 新型算法的特性分析

通過前面兩個(gè)算例的計(jì)算分析,驗(yàn)證了本文提出的降階模型新型算法是準(zhǔn)確有效的,特別是在計(jì)算效率方面,較常規(guī)算法有大幅度的提升.

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

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

4 結(jié)論

本文針對導(dǎo)熱系數(shù)隨溫度變化的一類瞬態(tài)非線性熱傳導(dǎo)問題,發(fā)展并提出了一種基于POD 模型降階技術(shù)的新型高效快速分析方法.通過單元預(yù)轉(zhuǎn)換方法和多級線性化方法的有效結(jié)合解決了非線性降階模型的常規(guī)計(jì)算方法在系數(shù)矩陣更新和迭代等計(jì)算環(huán)節(jié)的費(fèi)時(shí)問題,從而能夠獲得很高的計(jì)算效率.

通過二維和三維算例分析,演示了瞬態(tài)非線性熱傳導(dǎo)降階模型新型計(jì)算方法的應(yīng)用細(xì)節(jié),通過與全階模型有限元法計(jì)算結(jié)果的比較,驗(yàn)證了該新型算法的準(zhǔn)確性,并通過與降階模型常規(guī)計(jì)算方法以及近年來發(fā)展出的其他新方法在計(jì)算時(shí)間上的比較,驗(yàn)證了該新型算法的高效性.

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

此外,采用常數(shù)熱邊界條件下得到的POD 模態(tài)建立的瞬態(tài)非線性熱傳導(dǎo)降階模型可以對相同計(jì)算域在各種復(fù)雜的光滑/非光滑時(shí)變熱邊界條件下的非線性熱傳導(dǎo)過程進(jìn)行快速、準(zhǔn)確的分析和預(yù)測,這進(jìn)一步拓展了本文所提降階模型新型算法在工程應(yīng)用中的價(jià)值.