一種求解瞬態(tài)熱傳導方程的無條件穩(wěn)定方法1)

季 奕 邢譽峰

* (北京航空航天大學固體力學研究所,北京 100083)

? (北京航空航天大學高等理工學院,北京 100083)

引言

時間積分方法被廣泛應用在瞬態(tài)響應的求解中.根據(jù)遞推公式的格式,其可被劃分為兩類,分別是顯式方法和隱式方法.比較而言,顯式方法效率高但條件穩(wěn)定,隱式方法效率低但多為無條件穩(wěn)定.對于非線性系統(tǒng),確定顯式方法的時間步長具有挑戰(zhàn)性.因此,本文致力于發(fā)展對線性和非線性系統(tǒng)都是無條件穩(wěn)定的隱式方法.

對瞬態(tài)熱傳導問題,廣義梯形法則[1]是一種較為流行的方法,由它可衍生出Crank-Nicolson 方法、Galerkin 方法和后向差分方法(backward difference formula,BDF) 等.在這些方法中,只有Crank-Nicolson 方法具有二階精度,但其在高頻段無法提供數(shù)值耗散,從而會誘發(fā)出虛假的數(shù)值振蕩.一階精確的BDF 具有極強的數(shù)值耗散從而能迅速地過濾掉高頻信息,但其無法準確地保留重要的低頻信息.相比之下,Galerkin 方法能夠有效地平衡低頻精度和高頻耗散.隨著人們對精度、效率、耗散和穩(wěn)定性的追求,在過去的幾十年中,不斷涌現(xiàn)出新的時間積分方法[2-31].

除減小步長外,提高算法階次[2-6]也是一種提高精度的有效手段.針對一階線性常微分方程,Fung借助加權殘量法[2]、最小二乘法[3]和微分求積法[4]構造出一系列具有任意高階精度的無條件穩(wěn)定時間積分方法,但這些方法在求解過程中需要對系統(tǒng)矩陣進行擴維處理,從而大幅增加計算量.為了同時獲得高精度和高效率,一種精細積分方法[7]被提出.該方法利用Taylor 級數(shù)來近似解析的傳遞矩陣,利用分段線性技術處理外載荷.此外,該方法采用了2m算法和增量矩陣存貯技術,大幅提高了效率、減少了舍入誤差.雖然精細積分方法是條件穩(wěn)定的,但由于在計算中采用了極小步長,因此條件穩(wěn)定性并不影響其應用.精細時間積分方法已被廣泛地應用在熱傳導[8-9]、雙邊值[10-11]和結構動力學[12-13]等問題中.

上面提到的優(yōu)秀方法僅能求解線性問題,下面回顧一下能夠求解一般非線性瞬態(tài)熱傳導問題的方法[14-31].Runge-Kutta 方法[14-16]是具有無條件穩(wěn)定的高階格式.但與單步方法相比,多級結構使其計算量更高.為了在不犧牲計算效率的前提下提高精度,線性多步方法[17-21]和多分步方法[22-28]被提出.一般來說,線性多步方法是二階精確的,但其不能自啟動.與單步方法相比,線性多步方法的使用略顯繁瑣.多分步方法是自啟動的,但各個分步的有效剛度矩陣可能不同,致使其計算量略高于單步方法.瞬態(tài)熱傳導問題普遍存在于航空航天、土木和冶金等領域中,構造一種具有無條件穩(wěn)定性且能夠靈活處理該類問題的高精度、高效率時間積分方法是必要的.在這種背景下,本文提出了一種能處理一般瞬態(tài)熱傳導問題的無條件穩(wěn)定單步時間積分方法,其具備二階精度、高頻耗散和自啟動特性.利用拉格朗日插值函數(shù)分別近似溫度場及其一階導數(shù)場,并利用加權殘量法建立二者之間的關系.理論分析和數(shù)值測試均顯示出了本文方法在精度、耗散和穩(wěn)定性方面的優(yōu)勢.

1 算法描述

瞬態(tài)熱傳導問題的非線性控制方程[32]可以寫為

式中,C和K分別是N×N的熱容矩陣和熱傳導矩陣,T和Q分別是N×1 的溫度列向量和外部施加的熱源列向量.當C和K是與溫度T無關的定常矩陣時,式(1)退化為線性方程,即

為求解式(1)和式(2)中給出的瞬態(tài)熱傳導問題,本文首先把待求的時間域[0,tTotal]等間距地劃分為n個時間單元,各離散點為0=t0

式中,ψj(t)是與第j個節(jié)點關聯(lián)的拉格朗日插值函數(shù),Tj和分別描述了與第j個節(jié)點關聯(lián)的溫度及其一階導數(shù).拉格朗日插值函數(shù)的定義為

其中,τ=(t?tk)/Δt,τj用來確定第j個節(jié)點的位置即tj=tk+τjΔt.在本文方法中,τ0=0,τ1=0.5,τ2=1.在式(3)中,利用拉格朗日函數(shù)作為基函數(shù)分別獨立得到了T和,因此二者之間不存在關系T˙ =dT/dt.下面利用加權殘量法來建立T和T˙ 之間的關系.定義如下加權殘量r(t)

對于解析方法,上述殘量在任意時刻均為0.為了最小化r(t),本文方法要求

式中,w(t)為權函數(shù).將式(5)代入式(6)可得

選用不同的權函數(shù)w(t),可以獲得不同數(shù)值特性的時間積分方法.但這種方式很難設計出具有無條件穩(wěn)定和理想耗散的方法,為此本文沒有直接使用權函數(shù)w(t),而是利用如下參數(shù)θk來控制算法的性能.

將式(4)和式(8)代入式(7)整理可得

本文方法的數(shù)值性能僅由參數(shù)θ1和θ2控制,在之后的章節(jié)中,它們將按所希望的性能要求進行設計.

2 改進Hughes 理論

經典Hughes[33]理論可用來評估已有時間積分方法在非線性瞬態(tài)熱傳導系統(tǒng)中的穩(wěn)定性特性.為了能夠設計出對該類非線性問題無條件穩(wěn)定的新方法,本文對經典Hughes 理論進行了改進.

但大量的數(shù)值實驗表明對線性系統(tǒng)無條件穩(wěn)定的方法對非線性系統(tǒng)可能是失效的.在用時間積分方法求解非線性方程(1)時,由于不同時刻的C和K是不同的,因此利用線性化方法得到的不同時刻的熱特征值和特征向量是不同的.為了能夠考慮非線性系統(tǒng)的這個特征,Hughes[33]建立了如下模型

其中 λtk+1為tk+1=tk+Δt=(k+1)Δt時刻的熱特征值.若系統(tǒng)是線性的,則 λtk+1不隨時間變化,恒為常數(shù),即 λtk+1=λ,此時模型(14)退化為經典線性模型(13).將時間積分方法的遞推公式代入模型(14)可推出時變的傳遞因子A.若|A|≤ 1,則時間積分方法對單自由度非線性系統(tǒng)是無條件穩(wěn)定的.對于非線性多自由度系統(tǒng),對其控制方程(1)在不同時刻分別解耦時,同一階次的φ在不同時刻是不同的.不同時刻的溫度坐標和特征向量的順序是根據(jù)對應時刻特征值λ的從小向大的順序依次排列的.值得指出的是,由模型(14)推出的穩(wěn)定性條件對單自由度非線性系統(tǒng)是充分且必要的,而對多自由度非線性系統(tǒng)是必要但不充分的.借助方程(14),Hughes[33]詳細分析了廣義梯形法則在非線性瞬態(tài)熱傳導系統(tǒng)中的穩(wěn)定性.廣義梯形法則[32]的遞推公式為

將式(15)代入式(14)中可得傳遞因子為

將式(16)代入|A|≤ 1 中可得

從式(17)中可以發(fā)現(xiàn),在廣義梯形法則中只有BDF (α=1)對非線性系統(tǒng)是無條件穩(wěn)定的,而對線性系統(tǒng)無條件穩(wěn)定的Crank-Nicolson 方法(α=1/2)對非線性系統(tǒng)則是條件穩(wěn)定的.另外,注意到對于顯式方法(α=0),失穩(wěn)原因是時間步長尺寸Δt不夠小,而隱式方法(0 <α≤ 1)失穩(wěn)的原因與熱特征值λ演化規(guī)律相關.為了能夠借助Hughes 理論設計出對非線性瞬態(tài)熱傳導系統(tǒng)是無條件穩(wěn)定的隱式時間積分方法,本文在Hughes 理論中引入一個用來刻畫λ演化規(guī)律的函數(shù),其定義為

函數(shù)δt可以刻畫λ在一個時間單元內[tk,tk+Δt]的變化情況.將函數(shù)δt代入模型(14)可得改進的模型為

利用改進的模型(19)得到的廣義梯形法則的穩(wěn)定條件為

比較式(20)和式(17)可以發(fā)現(xiàn),由改進模型得到的穩(wěn)定性條件與Hughes 理論得到的結果完全一致,從而說明改進理論在穩(wěn)定性分析方面與Hughes理論是等價的.下面以廣義梯形法則的穩(wěn)定性條件(20)為例,簡述如何利用改進Hughes 理論設計無條件穩(wěn)定方法.從式(20)可以看出,當(? δtk+1α-α+1)≤0 成立時,不等式(20)恒成立.為了消除正的 δtk+1對穩(wěn)定性的影響,α的取值區(qū)間應為α≥1.當α≥1 時,廣義梯形法則是無條件穩(wěn)定的,其中α=1 對應的是BDF.這種通過消除 δtk+1對穩(wěn)定性影響的方法可以用來設計無條件穩(wěn)定方法,并被用于本文方法的參數(shù)設計中.Hughes 理論原模型(14)雖然可以用于分析方法的穩(wěn)定性,但難以用于設計無條件穩(wěn)定方法.

3 算法設計與分析

將本文方法的遞推公式(9)代入式(19)可得遞推關系為

其中,時變的傳遞因子A和載荷算子L分別為

式中,Ω= λtkΔt,函數(shù)gi,和hi(i=1,2)為

其中 δtk+1/2= λtk+1/2/ λtk,δtk+1= λtk+1/ λtk,tk+1/2=tk+Δt/2=(k+1/2)Δt.

首先考慮本文方法的高頻耗散特性.若一個無條件穩(wěn)定時間積分方法的傳遞因子滿足A|Ω→∞=0,則其被認為是L 型耗散方法.為了能夠快速過濾掉高頻振蕩,L 型耗散是被期望的,為此要求

從而可以建立θ1和θ2的約束關系為

進而式(22)中的傳遞因子A和式(23)中的載荷算子L分別被更新為

下面討論本文方法的穩(wěn)定性.這里假設1/2≤θ1≤3/4,從而使得式(27)中的傳遞因子A的分母恒為正值,進而不等式|A|≤ 1 可等價表示為

式中,Ω∈[0,+∞),δtk+1/2∈(0,+∞),δtk+1∈(0,+∞).從式(29)中可以觀察到,當0 < δtk+1<1 時,算法有可能失穩(wěn).為了消除這個不穩(wěn)定因素,令θ1=3/4.另外可以注意到,對于線性系統(tǒng)(δtk+1/2≡ δtk+1≡ 1),不等式(29)是恒成立的.再次說明造成隱式方法對非線性問題失穩(wěn)的原因與熱特征值的時變性有關.至此,兩個自由參數(shù)θ1和θ2確定完成,它們使得本文方法具有無條件穩(wěn)定性和L 型耗散.

最后,利用局部截斷誤差σ[32]來分析本文方法的精度階次.局部誤差定義為

將式(27)和式(28)代入上式整理可得

可以看到對線性系統(tǒng)(δtk+1/2≡ δtk+1≡ 1),s0=s1=0,故本文方法具有嚴格的二階精度.對于非線性系統(tǒng),當 δtk+1/2和 δtk+1均趨近于1 時,s0和s1則趨于0,此時本文方法是近似二階精確的.

為便于讀者使用,下面給出本文方法對線性系統(tǒng)(2)的求解流程.求解過程中涉及tk+1/2和tk+1兩個時刻的平衡方程,即

首先將遞推公式(9)代入平衡方程(33)中可解出tk+1和tk+1/2兩個時刻的溫度,即

將求得的Ttk+1和Ttk+1/2回代到式(9) 即可得到,從而獲得所有未知的變量信息.計算步驟如下:

第1 步準備工作:

(1) 構造熱容矩陣C,熱傳導矩陣K,熱源向量Q.

(2) 初始化T0和T˙0.

(3) 選擇步長Δt和計算算法參數(shù):

(4) 構造系數(shù)矩陣:

第2 步迭代運算:

(1) 計算tk+Δt時刻的有效載荷向量:

(2) 求解tk+1和tk+1/2時刻的溫度:

(3) 求解tk+Δt時刻溫度的一次導數(shù):

對非線性系統(tǒng)(1),利用遞推公式(9)、平衡方程(33)和Newton 迭代法,可求解出未知變量Ttk+1,Ttk+1/2和.

4 數(shù)值性能

本節(jié)對本文方法的無條件穩(wěn)定性、L 型耗散和二階精度進行數(shù)值驗證和討論.

4.1 傳遞因子

圖1(a)給出了本文方法和一些現(xiàn)有方法在線性系統(tǒng)中的傳遞因子[32]曲線,其中解析傳遞因子的表達式為Aexact=exp(?λΔt).與Crank-Nicolson 方法和Galerkin 方法相比,本文方法具有高頻耗散優(yōu)勢.與BDF 相比,本文方法有低頻精度優(yōu)勢.圖1(b)給出了本文方法在非線性系統(tǒng)下的傳遞因子曲線,可以看到對任意 δtk+1/2和 δtk+1的組合,在不同Ω值下|A|≤1 始終成立,并且?guī)缀鯚o數(shù)值振蕩區(qū).隨著Ω值增加,A逐漸趨于0,說明本文方法對非線性問題具備較強的高頻耗散能力.

圖1 傳遞因子Fig.1 Amplification factor

4.2 收斂率

收斂率可檢測時間積分方法的精度階次,其定義為在指定時刻物理量的相對誤差隨步長減小而減小的速率.先測試線性系統(tǒng)[34],考慮如下方程

圖2(a)給出了t=100 時的相對誤差,可以看出本文方法與Crank-Nicolson 方法具有相同的斜率.從而說明對線性系統(tǒng),本文方法具有嚴格的二階精度.再測試非線性系統(tǒng)[34],考慮如下方程

該問題的解析解為T(t)=T0/(3σT03t+1)1/3.圖2(b)給出了t=10 的相對誤差,可以看出對非線性問題,本文方法也能夠達到二階精度,此時s0和s1趨于0.

圖2 收斂率Fig.2 Convergence rate

5 數(shù)值測試

本節(jié)利用3 個算例來測試本文方法在不同類型瞬態(tài)熱傳導問題中的性能.本文方法的求解過程中涉及tk+1和tk+1/2兩個時刻的平衡方程,而廣義梯形法則僅求解tk+1時刻平衡方程.為保證精度比較在相同計算量下執(zhí)行,在所有的算例中,本文方法的時間步長取為單步方法的兩倍.在本節(jié),除廣義梯形法則外,具有二階精確和L 型耗散的BDF2 方法[21]也被考慮,由于其不具備自啟動特性,本文采用Crank-Nicolson 方法作為啟動算法.此外,所有問題的參考解均由使用小步長的Crank-Nicolson 方法得到.

5.1 矩形功能梯度板中的二維熱傳導問題

考慮一個功能梯度板中的二維熱傳導問題[35],如圖3 所知,假設導熱系數(shù)kx1=kx2=1+(x1? 1)2+(x2? 1)2,密度ρ=1,比熱容cv=1.整個域內施加單位初始溫度,所有邊界施加T=0 的溫度場.該矩形域利用20×20 個4 結點矩形單元進行空間離散,廣義梯形法則采用步長Δt=0.01.

圖3 算例1 的幾何尺寸和邊界條件Fig.3 Geometry and boundary conditions for Ex.1

圖4 給出了中點A(1,1)和角點B(1.9,1.9)的溫度?時間曲線.可以看到在中點A處,所有方法均與參考解吻合得較好,但在邊界點B處,Crank-Nicolson 方法存在明顯的數(shù)值振蕩.

圖4 點A 和點B 的溫度?時間曲線Fig.4 Temperatures of point A and point B versus time

為了進一步比較這些方法在域內的精度表現(xiàn),圖5 提供了所有方法在t=0.1 時刻溫度誤差的分布結果,可以看出本文方法不會像Crank-Nicolson 方法一樣在邊界處產生劇烈的數(shù)值振蕩,且在所有方法中,本文方法的域內精度最高.BDF 和Galerkin 方法雖然沒有數(shù)值振蕩,但域內精度不令人滿意.

圖5 溫度誤差分布(t=0.1)Fig.5 Errors of temperature at t=0.1

5.2 一維非線性熱傳導問題

考慮一個材料參數(shù)與溫度有關的一維桿[36],桿長l假設為1,材料參數(shù)為k=γT2+βT+η(γ,β,η∈R)和ρcv=1.在桿的兩端施加T=0 的溫度場,域內施加單位初始溫度.利用20 個兩結點單元對桿進行空間離散.

首先考慮γ=0,β=η=1 的情況,廣義梯形法則的Δt=0.000 25.圖6 給出了中點處(x=0.5)和邊界處(x=0.95)的溫度?時間曲線,可以看到本文提出的方法展示出明顯的精度優(yōu)勢.

圖6 溫度?時間曲線(γ=0,β=η=1)Fig.6 Temperature-time history

圖6 溫度?時間曲線(γ=0,β=η=1) (續(xù))Fig.6 Temperature-time history (continued)

為了展示本文方法的穩(wěn)定性優(yōu)勢,另一種情況β=100,γ=η=1 被考慮.此時廣義梯形法則的Δt=0.005.從圖7 可以看到Crank-Nicolson 方法在計算一段時間后就失去了穩(wěn)定性,而本文方法不僅沒有失去穩(wěn)定性,且在所有收斂的方法中具有最高的精度.

圖7 溫度?時間曲線(β=100,γ=η=1)Fig.7 Temperature-time history

5.3 具有外部熱源的二維熱傳導問題

考慮一均質矩形板[9],如圖8 所示,外部熱源f(t)=cos(10t)從右邊界持續(xù)不斷地向域內輸入,底部邊界施加T=0 的溫度場,其余邊界為絕熱狀態(tài).利用50×50 個4 結點矩形單元進行劃分,廣義梯形法則的Δt=0.005.

圖8 算例3 的幾何尺寸和邊界條件Fig.8 Geometry and boundary conditions for Ex.3

圖9 給出了點A(0.98,0.98)處的溫度和溫度一階導數(shù)的結果,從絕對誤差的數(shù)值上可以看出:本文方法要比二階Crank-Nicolson 方法和BDF2 更加精確.總體來說,這3 種二階精度算法的精度明顯優(yōu)于一階精度的Galerkin 方法和BDF.

圖9 點A 處溫度和其一階導數(shù)與時間的關系曲線Fig.9 Temperature and its derivative -time histories of point A

6 結論

本文提出了一種適用于求解一般瞬態(tài)熱傳導問題的無條件穩(wěn)定單步時間積分方法.該方法具有無條件穩(wěn)定性、二階精度、L 型數(shù)值耗散和自啟動特性.需要強調的是,與大多數(shù)現(xiàn)有時間積分方法不同,本文方法對線性和非線性瞬態(tài)熱傳導問題均是無條件穩(wěn)定的.

線性和非線性數(shù)值比較結果表明,與著名的二階Crank-Nicolson 方法相比,在計算量相同的前提下,本文方法具有更理想的精度和穩(wěn)定性,并且不會在高頻段誘發(fā)出虛假的數(shù)值振蕩.與具有L 型耗散的二階BDF2 相比,本文方法具有更高的精度且無需其他算法來啟動計算.

鑒于本文方法具有優(yōu)秀的數(shù)值性能和易執(zhí)行性,故推薦其用于一般瞬態(tài)熱傳導問題的求解.