基于特征分裂有限元準(zhǔn)隱格式的共軛傳熱整體耦合數(shù)值模擬方法1)

劉 瑜 鄧家鈺 王成恩,2) 蘇紅星

?(上海交通大學(xué)機(jī)械與動力工程學(xué)院,上海 200240)

?(大連理工大學(xué)航空航天學(xué)院,遼寧大連 116024)

引言

共軛傳熱現(xiàn)象在科學(xué)和工程領(lǐng)域中大量存在[1],比如核反應(yīng)系統(tǒng)[2],燃?xì)廨啓C(jī)渦輪葉片冷卻[3-4],航天器熱防護(hù)[5-6],化工處理[7],生物醫(yī)學(xué)[8-9],電子設(shè)備散熱[10],MEMS 設(shè)備中的微槽道對流換熱[11]等.

共軛傳熱問題除了實驗研究以外,可以通過理論分析或數(shù)值模擬的方法求解.理論分析方法對簡單的模型(如平板,圓管等)進(jìn)行建模,通過簡化得到精確或近似的理論解[12].理論方法對共軛傳熱問題的理解很重要,但只適用于簡單的問題.數(shù)值模擬將流動的數(shù)值解與固體傳熱的數(shù)值解耦合起來,主要有兩種耦合方法:分區(qū)耦合(partitioned)[13]和整體耦合(monolithic)[14]方法.分區(qū)耦合方法對流場和固體熱分析采用各自相應(yīng)的成熟的求解器,然后通過迭代方法滿足流熱耦合的界面條件,即在界面上溫度和熱流應(yīng)是連續(xù)的.分區(qū)耦合方法的優(yōu)點是可以采用現(xiàn)有的流動和固體熱分析求解程序,但是為了滿足界面條件需要進(jìn)行迭代求解,在不同的求解器之間傳遞邊界條件,算法實現(xiàn)復(fù)雜,并且還需要考慮分區(qū)耦合算法在流固耦合界面上的穩(wěn)定性條件[15-20].整體耦合方法將流動方程和固體傳熱方程統(tǒng)一離散,然后求解單一的方程.整體耦合方法的優(yōu)點是不需要處理流固耦合界面條件,缺點是如果流動特征時間和固體傳熱特征時間差別很大,非定常共軛傳熱的計算量將是巨大的[21-22].

目前存在大量的數(shù)值方法模擬共軛傳熱問題，如有限差分法[23],有限體積法[6,24-25],浸入式邊界法[26],無網(wǎng)格法[27],格子玻爾茲曼方法[28-30],邊界元法[31],有限元法等.有限元法在計算固體傳熱問題上具有優(yōu)勢，如果流場求解也采用有限元求解,可以降低程序?qū)崿F(xiàn)的復(fù)雜性.采用有限元法求解共軛傳熱的文獻(xiàn)相對較少.Misa 和Sarkar[32]采用整體耦合方法,利用滿足inf-sup 條件的經(jīng)典伽遼金(Galerkin)有限元方法對方腔自然對流的共軛傳熱進(jìn)行了模擬.Malatip 等[33]采用Rice 和Schnipke[34]提出的交錯有限元方法和整體耦合方法求解定常共軛傳熱問題.Malatip 等將Choi 和Yoo[35]提出的二階時間精度分裂有限元方法推廣到求解共軛傳熱問題[36]和流?熱?結(jié)構(gòu)耦合問題[37].Verstraete 和Scholl[19]以及Scholl 等[20]采用壓力穩(wěn)定有限元法計算流動,定常有限元法計算固體傳熱,對分區(qū)耦合方法計算定常共軛傳熱的穩(wěn)定性問題進(jìn)行了研究.Long 等[38]將光滑有限元方法(smoothed finite element method) 和改進(jìn)的光滑粒子流體動力學(xué)方法耦合起來求解熱?流?結(jié)構(gòu)相互作用問題.

利用伽遼金有限元法求解不可壓縮流動,需要有限單元滿足所謂的inf-sup 條件,即速度和壓力的插值函數(shù)不能具有相同的階數(shù)[39].這增加了算法實現(xiàn)的難度,特別是對于三維問題更是如此.為了采用等階元計算,Hughes 等[40]提出了用PSPG (Pressure-Stabilizing/Petrov-Galerkin) 穩(wěn)定項來規(guī)避inf-sup 條件,為了實現(xiàn)對流主導(dǎo)流動的穩(wěn)定計算還需要添加SUPG(Streamline-Upwind/Petrov-Galerkin) 穩(wěn)定項[41].Zienkiewicz 等[42]在分?jǐn)?shù)步法的基礎(chǔ)上,提出了基于特征分裂的有限元法(CBS).CBS 方法有半隱格式(semi-implicit),準(zhǔn)隱格式(quasi-implicit) 和基于人工壓縮性的格式(CBS-AC)[42].準(zhǔn)隱格式的時間步長只取決于對流穩(wěn)定性,計算效率高于半隱格式和基于人工壓縮性的格式[43].在作者閱讀的有限文獻(xiàn)范圍內(nèi),還沒有看到將CBS 方法的準(zhǔn)隱格式用于共軛傳熱數(shù)值模擬的例子.CBS 方法既適用于不可壓縮流動的求解,也適用于亞聲速、超聲速流動的求解;并且同樣可以采用等階插值函數(shù);CBS 方法是分?jǐn)?shù)步方法,每一個分裂步得到的方程是標(biāo)量方程,對存儲量要求小,方程的求解也相對容易.本文的主要目的是采用CBS 方法的準(zhǔn)隱格式,發(fā)展一種整體耦合層流共軛傳熱數(shù)值模擬方法.

本文接下來給出共軛傳熱的控制方程,并介紹流場計算采用的基于特征分裂的有限元方法,隨后給出整體耦合共軛傳熱數(shù)值模擬算法.為了驗證本文所采用算法的準(zhǔn)確性,對不可壓縮流動問題和共軛傳熱問題進(jìn)行了模擬.

1 計算方法

1.1 控制方程

非守恒形式的不可壓縮Navier-Stokes 方程為

上式是無量綱化的方程,其中ui,(i=1,2,3)是流體速度分量,p是流體靜壓,Re為雷諾數(shù).方程中變量下標(biāo)滿足愛因斯坦求和規(guī)則.

假設(shè)流體的密度是常量,且其他流體性質(zhì)也與溫度無關(guān).將能量方程表示為溫度的形式,即

其中,T為溫度,Qf為流體中的熱源項,Pr為普朗特數(shù).在這種情形下,能量方程與連續(xù)性方程和動量方程解耦.溫度只是簡單的輸運變量,在得到了速度場后,計算溫度方程就可以得到溫度場的分布.

固體傳熱方程為

其中κsf=κs/κf,cr=ρscps/ρfcpf.ρf,ρs分別是流體和固體的密度,κf,κs分別是流體和固體的導(dǎo)熱系數(shù),cps和cpf分別表示固體和流體的比熱.推導(dǎo)以上無量綱方程采用的參考量見附錄.

1.2 求解不可壓縮流動的分?jǐn)?shù)步法

方程(1)的半離散形式為

其中?t是時間步長,θ1,θ2和θ3分別是控制對流項、壓力項和黏性項時間離散的參數(shù),其取值范圍為[0,1].為了避免求解非線性方程,只考慮θ1=0 的情形.在推導(dǎo)CBS 格式之前,首先引入輔助變量?u?和?u??,且令

采用基于Helmholtz-Hodge 分解的2 級分?jǐn)?shù)步方法,將方程分裂為

1.3 對流穩(wěn)定項

對于對流占優(yōu)問題,CBS 格式需要添加對流穩(wěn)定項[43].本文在半離散方程(5) 中添加顯式的穩(wěn)定項,即

其中?ts為對流穩(wěn)定項需要的時間步長,可以認(rèn)為是穩(wěn)定性參數(shù),與SUPG 中的穩(wěn)定參數(shù)τ[41]的作用相似.在分?jǐn)?shù)步方法中,穩(wěn)定項也要進(jìn)行分裂,在預(yù)測步式(7)中添加的穩(wěn)定項為

即不考慮壓力項.在修正步式(8) 中,穩(wěn)定項僅包含壓力項,即

類似可以得到帶穩(wěn)定項的能量方程半離散形式,即

其中最后一項是對流穩(wěn)定項.

1.4 基于特征分裂有限元方法的矩陣形式

對半離散方程采用伽遼金有限元方法進(jìn)行空間離散,根據(jù)參數(shù)的不同,可以得到不同的格式[43].

1.4.1 半隱格式

步驟1

其中Mu,i,C,K分別表示速度質(zhì)量矩陣,離散對流算子和黏性算子.Fi表示黏性邊界積分,fs,i是穩(wěn)定項.

步驟2

步驟3

其中L,D,G分別表示離散拉普拉斯算子,離散散度算子和離散梯度算子,fp,i是穩(wěn)定項.半隱格式是條件穩(wěn)定的,穩(wěn)定性由對流項和黏性項確定.

在計算得到了速度場后,即可以將速度場代入到能量方程計算溫度場.

步驟4

其中,MT,KT分別是溫度方程質(zhì)量矩陣和溫度擴(kuò)散算子,FT表示溫度擴(kuò)散項邊界積分,QT為熱源項,fT為穩(wěn)定項.

半隱格式的單元穩(wěn)定時間步長為

其中上標(biāo)e表示單元,h是單元的網(wǎng)格尺度,?tcon,?tvis分別是對流時間步長和黏性時間步長.全局時間步長取所有單元穩(wěn)定時間步長的最小值,即

黏性穩(wěn)定性在一定條件下會對時間步長產(chǎn)生嚴(yán)重的限制.為了克服這一缺點,提出了準(zhǔn)隱格式.

1.4.2 準(zhǔn)隱格式

與半隱格式不同,準(zhǔn)隱格式對黏性項進(jìn)行隱式處理[43],即1/2θ31.

步驟1

黏性項的邊界積分仍然作顯式處理.步驟2 和步驟3與半隱格式相同.

步驟3 還可以采用非投影形式,即對于速度修正步

步驟4

準(zhǔn)隱格式的時間步長只受對流項限制,即

以上矩陣和源項的具體形式見附錄.

說明1:半隱格式的時間步長受黏性穩(wěn)定性條件影響;而基于人工可壓縮性的格式(CBS-AC),根據(jù)文獻(xiàn)[43]報道,準(zhǔn)隱格式的計算效率比其高3 ～100 倍.因此本文計算采用準(zhǔn)隱格式.

說明2:CBS 方法允許等階插值.本文選擇等階有限元,即速度、壓力、溫度都采用相同的有限元基函數(shù),這大大方便了程序的實現(xiàn).

說明3:準(zhǔn)隱格式求解的四步方程的矩陣是對稱正定的,故可以采用共軛梯度法求解.預(yù)處理矩陣本文采用雅克比方法計算.

說明4:本文計算所采用的參數(shù)為θ1=0,θ2=θ3=θ4=1.

1.5 共軛傳熱耦合算法

在進(jìn)行流固耦合模擬時,需要劃分流場網(wǎng)格和固體網(wǎng)格.本文設(shè)計了一種流固耦合網(wǎng)格和邊界條件表達(dá)方式.出于簡單考慮,流固耦合界面是適配的,即在流固耦合界面上流場網(wǎng)格點和固體網(wǎng)格點是一致的.為了方便,流體區(qū)域和固體區(qū)域作為一個統(tǒng)一的計算域劃分單一的網(wǎng)格.為了區(qū)分流體和固體網(wǎng)格,需要用一個指標(biāo)集來表示網(wǎng)格單元的材料屬性,考慮用一個材料屬性文件來存儲網(wǎng)格單元的屬性,即mateiral.txt.

對于溫度場而言,流固耦合邊界屬于內(nèi)面,因此在流固耦合邊界上只需要指定流場變量邊界條件.而在固體的其他邊界上只需要指定溫度邊界條件,不需要指定流場變量的邊界條件.

如果采用前面提及的流固耦合網(wǎng)格及邊界條件,可以得到共軛傳熱耦合算法如Algorithm 1 所示.

固體傳熱方程離散后的形式為

溫度方程在流體區(qū)域按照式(27) 離散,在固體區(qū)域按照式(29) 離散,最后將溫度方程裝配為統(tǒng)一的方程組.

該算法是一種整體耦合算法,與分區(qū)算法相比,不需要通過迭代就可以隱式滿足流熱耦合界面條件.

1.6 時間步長

1.6.1 流體時間步長

Bevan 等[43]研究了CBS 方法中的半隱格式、準(zhǔn)隱格式和基于人工壓縮性的格式的計算效率,發(fā)現(xiàn)準(zhǔn)隱格式由于只需要考慮對流時間步長,相比于其他兩種格式具有較高的計算效率.文獻(xiàn)[43] 中沒有對時間推進(jìn)步長和穩(wěn)定項中的時間步長進(jìn)行區(qū)分,在穩(wěn)定項中采用的也是時間推進(jìn)步長,即?ts=?t,并且建議采用半隱格式和準(zhǔn)隱格式時,時間步長應(yīng)該選擇全局時間步長.

作者認(rèn)為穩(wěn)定項中的時間步長?ts可以視為穩(wěn)定性參數(shù),而穩(wěn)定性參數(shù)應(yīng)具有當(dāng)?shù)匦再|(zhì).當(dāng)計算采用全局時間步長時,如果整個計算區(qū)域的穩(wěn)定性時間步長都取同樣的值,有可能不能較好地抑制數(shù)值不穩(wěn)定性.因此本文采用準(zhǔn)隱方法進(jìn)行計算,對物理推進(jìn)時間步長和穩(wěn)定項中的時間步長做了區(qū)分.物理推進(jìn)時間步長?t按式(28)計算，取全場的最小時間步長作為全局時間步長.穩(wěn)定項中的時間步長為局部時間步長?ts，按照式(22)計算.

1.6.2 固體時間步長

在固體傳熱方程(4)中,固體密度與比熱的乘積與流體密度與比熱的乘積之比,cr,在一般情況下有cr?1.這導(dǎo)致傳熱方程時間尺度遠(yuǎn)大于流動方程的時間尺度.He 等[44]估計了典型共軛傳熱問題的固體和流體的穩(wěn)定時間步長之比為103～105.當(dāng)對共軛傳熱問題的流體溫度方程和固體溫度方程統(tǒng)一求解時,如果采用流體部分的穩(wěn)定時間步長進(jìn)行計算,則固體部分的溫度場需要迭代很多步才能收斂.對于定常問題,為了提高收斂速度,在固體部分需要采用比流體時間步長更大的時間步長

其中?tflu和?tsol分別表示流動計算和固體傳熱計算采用的時間步長,n為放大因子.

式(30)中如果選擇n=cr,則相當(dāng)于在方程(4)中令cr=1,?t=?tflu=?tsol.這是文獻(xiàn)[23,45]采用的做法.本文定義名義上的cr值為cr′,即

計算時用cr′代替cr.名義上固體時間步長與流體時間步長相同,但是通過調(diào)整cr′實際上改變了固體時間步長.cr′可以足夠小,在滿足穩(wěn)定性的前提下提高收斂速度.這將在下面的算例中得到驗證.

2 算法驗證

2.1 穩(wěn)定項時間步長選取方法對計算結(jié)果的影響

首先對穩(wěn)定項中的時間步長的選擇方式對計算結(jié)果的影響進(jìn)行比較.采用3 種方式計算穩(wěn)定項中的時間步長,分別是:與全局時間步長相同,采用局部時間步長和取零值(即不添加穩(wěn)定項).計算采用的例子是方腔頂蓋驅(qū)動流.方腔的頂蓋按指定的速度運動,方腔的其他邊界固定.在頂蓋的作用下,方腔中流體將發(fā)生運動.計算域? 為(0,1)×(0,1),邊界條件都為Dirichlet 條件.頂蓋的運動速度,按照文獻(xiàn)[46],給定如下

當(dāng)網(wǎng)格足夠稠密時,穩(wěn)定項的作用不明顯,為此在這里特意選擇比較稀疏的網(wǎng)格進(jìn)行計算.網(wǎng)格包含800個線性三角形單元和441 個節(jié)點.網(wǎng)格如圖1(a) 所示.圖1(b)～圖1(d) 分別給出了由3 種方法計算的Re=5000 時x方向上的速度分量u的云圖.圖2 給出了在y=0.8 上u的分布.可見如果不加入穩(wěn)定項,流場變量將出現(xiàn)很大的偽振蕩;采用文獻(xiàn)[43] 中的方法,由于穩(wěn)定項中的時間步長全場都是相同的,不能完全抑制偽振蕩的出現(xiàn);而采用本文的方法流場變量分布光滑,抑制了偽振蕩的出現(xiàn).

圖1 頂蓋驅(qū)動流的計算網(wǎng)格和速度云圖Fig.1 Computational mesh and contours of x-component of velocity

圖2 由不同穩(wěn)定項計算方法得到的速度的x 分量在y=0.8 上的分布Fig.2 Velocity profile at y=0.8 predicted by different stabilization schemes

2.2 層流計算驗證

2.2.1 方腔頂蓋驅(qū)動流

方腔頂蓋驅(qū)動流是用來檢驗層流計算準(zhǔn)確性的典型算例.計算和邊界條件與上節(jié)相同.采用一階三角形單元計算,單元數(shù)為3200,節(jié)點數(shù)為1681,在邊界附近進(jìn)行了加密.分別對Re=400,1000,5000,10 000 時的情形進(jìn)行了計算.圖3 比較了x=0.5 上速度分量v和y=0.5 上速度分量u的分布與Ghia等[47]的計算結(jié)果,可見符合很好.

2.2.2 層流后臺階流動

對于這個問題有實驗數(shù)據(jù)可供比較,因而為很多作者所引用.本文選擇雷諾數(shù)為229 的情形進(jìn)行計算.問題的計算域和邊界條件如圖4 所示.進(jìn)口速度u的分布的表達(dá)式為

在壁面處按無滑移邊界條件處理,在出口速度指定Neumann 條件,壓力指定為0,在其余邊界上壓力指定為Neumann 條件.計算網(wǎng)格包含15 200 個三角形單元,7841 個節(jié)點.

圖5 分別給出了壓力云圖和速度云圖.圖(6)給出了不同截面上速度分布與實驗數(shù)據(jù)[48]的比較,可見計算結(jié)果與實驗符合很好.

圖3 不同雷諾數(shù)條件下的速度型面Fig.3 Velocity profiles in the x and y directions

圖4 層流后臺階流動計算域和邊界條件Fig.4 Domain and boundary conditions of flow over a backward facing step

圖5 層流后臺階流動的壓力和速度云圖Fig.5 Pressure and velocity contours of flow over a backward facing step

2.3 共軛傳熱數(shù)值模擬

本文對若干共軛傳熱問題進(jìn)行模擬,以檢驗算法和程序的準(zhǔn)確性.其中,反向流動熱交換器在模擬中考慮了cr′的取值對收斂速度的影響,其余算例都設(shè)置cr′=1.

2.3.1 共軛傳熱庫埃特流

共軛傳熱庫埃特流[33]的示意圖為圖7.上壁以固定速度運動,溫度為T2=0;下壁固定具有厚度,需要考慮熱傳導(dǎo),下壁底部溫度為T1=1;流固耦合界面上溫度的精確解為

其中h1,h2分別為固體和流體區(qū)域的厚度.流體域中的溫度和固體域中的溫度在y方向上為線性分布.

圖8 給出了取不同的κsf值時,由一階三角形單元計算的溫度沿y方向上的分布曲線,網(wǎng)格單元數(shù)為1200.可見計算值與理論值十分吻合.

圖8 不同條件下共軛傳熱庫埃特流溫度分布曲線與理論解的比較Fig.8 Comparison of benchmark solutions for Couette flow condition

2.3.2 二維反向流動熱交換器

二維反向流動熱交換器的示意圖如圖9 所示.在熱交換器中分布兩個平行流道,流道被一金屬平板分隔開.兩個流道的厚度分別為δ1,δ3,隔板的厚度為δ2.流道的外壁是絕熱的,流道入口處的速度和溫度被認(rèn)為是均勻的.計算采用的參數(shù)如下:幾何參數(shù)δ1=δ2=δ3=0.1 m,L=1 m;流動參數(shù)U1=0.2 m/s,T1=800 K,U3=0.1 m/s,T3=300 K,ρf=1000 kg/m3,κf=10 W/(m·K),cpf=25 J/(kg·K),μ=0.15 kg(m·s);金屬板的物性ρs=8000 kg/m3,κs=50 W/(m·K),cps=500 J/(kg·K).

圖9 二維反向流動熱交換器Fig.9 A conjugate counter flow heat exchanger

計算網(wǎng)格采用線性三角形單元,單元數(shù)為3360,節(jié)點數(shù)為1763.

在本例中cr的準(zhǔn)確值為160.為了考慮cr′的取值對收斂速率的影響,分別對cr′采用準(zhǔn)確值和cr′=1,0.001 的情形進(jìn)行了模擬,圖10 給出了溫度殘值的收斂曲線.可見cr′如果采用準(zhǔn)確的值,收斂將非常緩慢.如果令cr′=1,收斂速度將大幅增大,繼續(xù)減小cr′的值，收斂速度會繼續(xù)增大.cr′的取值不影響最終的收斂解.

圖10 采用不同的cr′ 的收斂速率比較Fig.10 Comparison of convergence rate for different cr′ values

計算的速度場和溫度場分別如圖11(a)和圖11(b)所示.本文還計算了不同的固體和流體熱傳導(dǎo)系數(shù)之比條件下的共軛傳熱.在x=L/2 上溫度的分布曲線如圖12 所示,并與Malatip 等[37]計算的結(jié)果進(jìn)行比較,可見與文獻(xiàn)給出的結(jié)果符合良好.

圖11 二維反向流動熱交換器x 方向上的速度云圖和溫度云圖Fig.11 Predicted x-component of velocity and temperature contours for a counter flow heat exchanger

圖12 分別取κsf=0.1,1,5,10,二維反向流動熱交換器的溫度在x= L/2 上的分布與文獻(xiàn)[37]的比較Fig.12 Temperature distributions at x= L/2 for κsf=0.1,1,5,10 of a counter flow heat exchanger compared with results from Ref.[37]

2.3.3 加熱固體的強(qiáng)迫對流冷卻

第3 個測試?yán)邮橇鲃恿鹘?jīng)兩塊平行平板間含有3 個矩形加熱固體的強(qiáng)迫對流冷卻問題.該問題首次由Davalath 和Bayazitoglu[23]提出用來模擬集成電路元件的冷卻.該問題的計算域和邊界條件如圖13 所示,矩形塊尺寸相同,且等間距分布.流體以完全發(fā)展的速度型面從左側(cè)進(jìn)入,從右側(cè)流出.3 個固體塊中的生成的熱可以假設(shè)為具有均勻分布的熱源項,本文中假設(shè)熱源項Qs=8.計算時取Pr=0.71,κsf=10,分別計算了雷諾數(shù)為Re=100,500,1000 的情形.計算網(wǎng)格如圖14 所示,網(wǎng)格中含有個10 900 個線性三角形單元,5708 個網(wǎng)格節(jié)點.圖15 給出了雷諾數(shù)為100 時的壓力、速度和溫度云圖.圖16 給出了固體?流體界面上溫度的分布,并與文獻(xiàn)[37]的結(jié)果進(jìn)行了比較,可見本文計算結(jié)果與文獻(xiàn)符合很好.

圖13 加熱固體的強(qiáng)迫對流冷卻問題的計算域和邊界條件Fig.13 Domain and boundary conditions of forced convection cooling across rectangular blocks

圖14 加熱固體的強(qiáng)迫對流冷卻問題的計算網(wǎng)格Fig.14 Computational mesh of forced convection cooling across rectangular blocks

圖15 雷諾數(shù)為100 時,加熱固體的強(qiáng)迫對流冷卻的壓力,速度幅度和溫度云圖Fig.15 Pressure,velocity magnitude and temperature contours for Re=100

圖16 雷諾數(shù)為100,500,1000 時固體?流體界面上溫度的分布與文獻(xiàn)[37]的比較Fig.16 Comparison of wall temperature distributions along solid-fluid interface of rectangular blocks for Re=100,500,1000 with Ref.[37]

3 結(jié)論

本文通過將物理推進(jìn)時間步長與穩(wěn)定項中的時間步長進(jìn)行區(qū)別,改進(jìn)了基于特征分裂的有限元法的準(zhǔn)隱格式,使其具有了更好的穩(wěn)定性.在此基礎(chǔ)上,利用改進(jìn)了的準(zhǔn)隱格式和整體耦合方法實現(xiàn)了不可壓縮流體和固體間共軛傳熱的數(shù)值模擬并通過與理論、實驗和文獻(xiàn)中的數(shù)值解進(jìn)行比較,驗證了本文提出方法的正確性.還研究了固體區(qū)域時間步長對定常共軛傳熱問題數(shù)值模擬收斂性的影響,發(fā)現(xiàn)選擇合適的固體時間步長可以顯著提高收斂速率.

本文的一個不足之處是沒有對穩(wěn)定時間步長對基于特征分裂的有限元法準(zhǔn)隱格式的穩(wěn)定性的影響進(jìn)行理論分析,在未來的工作中將對此進(jìn)行研究;同時還可以研究物理推進(jìn)時間步長采用局部時間步長對定常共軛傳熱問題計算收斂性的影響.未來還將采用本文提出的方法對自然對流和混合對流問題進(jìn)行研究,進(jìn)一步探索該方法的適用范圍.

附錄

控制方程的無量綱化

有量綱形式的控制方程為: