溫度場(chǎng)的有限差分計(jì)算
——以一維溫度場(chǎng)為例①

趙銳，楊名揚(yáng)，淡丹輝，3

1.新疆大學(xué) 建筑工程學(xué)院， 烏魯木齊 830017； 2.新疆建筑結(jié)構(gòu)與抗震重點(diǎn)實(shí)驗(yàn)室， 烏魯木齊 830017； 3.同濟(jì)大學(xué) 橋梁工程系， 上海 200092

熱量的傳遞過(guò)程是通過(guò)熱傳導(dǎo)、熱對(duì)流和熱輻射3種方式進(jìn)行的， 其中， 熱傳導(dǎo)是通過(guò)固體或者靜止流體傳播熱量； 熱對(duì)流局限于液體和氣體， 并且往往涉及流體固體邊界的熱交換； 熱輻射能量是通過(guò)電磁波傳遞的[1]. 在計(jì)算物體溫度場(chǎng)時(shí)， 固態(tài)物體內(nèi)部熱量的傳遞方式為熱傳導(dǎo)， 將熱對(duì)流和熱輻射作為邊界條件. 分析傳熱問(wèn)題時(shí)， 從不同的角度出發(fā)可將傳熱問(wèn)題分為3類： ①根據(jù)熱量傳遞的維度可將熱傳導(dǎo)問(wèn)題分為一維熱傳導(dǎo)問(wèn)題、二維熱傳導(dǎo)問(wèn)題和三維熱傳導(dǎo)問(wèn)題； ②根據(jù)是否考慮溫度場(chǎng)隨時(shí)間變化， 可分為瞬態(tài)溫度場(chǎng)和穩(wěn)態(tài)溫度場(chǎng)； ③根據(jù)是否考慮材料性能隨溫度的變化， 可分為變系數(shù)熱傳導(dǎo)問(wèn)題和常系數(shù)熱傳導(dǎo)問(wèn)題.

熱量傳遞引起的溫度場(chǎng)變化涉及到生物[2]、醫(yī)療[3]、材料[4-6]， 建筑結(jié)構(gòu)[7-9]等多個(gè)領(lǐng)域， 在溫度場(chǎng)研究中， 將與溫度分布相關(guān)的物理參數(shù)稱為熱工參數(shù)， 包括比熱、熱傳導(dǎo)系數(shù)以及密度， 其中密度隨溫度變化微小. 對(duì)于特殊材料或環(huán)境溫度變化范圍較小的情況， 如自然環(huán)境等， 可不考慮材料熱工參數(shù)隨溫度的變化. 然而， 大多材料的熱工參數(shù)都會(huì)隨著溫度的變化而改變， 以鋼材為例， 各國(guó)規(guī)范均給出了鋼材熱工參數(shù)隨溫度的變化規(guī)律[10]. 為了求解熱工參數(shù)隨溫度變化的熱傳導(dǎo)問(wèn)題， 國(guó)內(nèi)外學(xué)者做了大量研究. 在穩(wěn)態(tài)溫度場(chǎng)的研究中， ünal[11]基于平面板模型提出一種考慮熱傳導(dǎo)率隨溫度線性變化的一維穩(wěn)態(tài)熱傳導(dǎo)問(wèn)題求解方法； Hussein[12]通過(guò)理論推導(dǎo)， 在球坐標(biāo)系中求解了實(shí)心球體在有內(nèi)部熱源、考慮材料熱傳導(dǎo)率隨溫度線性變化的情況下的一維穩(wěn)態(tài)傳熱問(wèn)題， 并且對(duì)比了使用3種不同材料時(shí)的溫度分布情況； 高仲儀[13]采用λ變換法求解變系數(shù)穩(wěn)態(tài)熱傳導(dǎo)問(wèn)題， 該方法假設(shè)熱傳導(dǎo)率隨溫度線性變化且熱傳導(dǎo)率與溫度一一對(duì)應(yīng)， 將溫度逆向表示為熱傳導(dǎo)率的函數(shù)， 通過(guò)求解熱傳導(dǎo)率進(jìn)而得到溫度場(chǎng)； 嚴(yán)治軍[14]給出了具有對(duì)流換熱邊界條件的一維常系數(shù)熱傳導(dǎo)差分解法， 考慮多層不同材料時(shí)溫度場(chǎng)分布， 最后采用FORTRAN語(yǔ)言編制成熱傳導(dǎo)解析的計(jì)算程序. 瞬態(tài)溫度場(chǎng)的研究也大多是基于一定假設(shè)簡(jiǎn)化得到的計(jì)算結(jié)果， 魏光坪[15]、康為江[16]和彭友松[17]在研究自然條件下的橋梁溫度場(chǎng)時(shí)， 由于自然環(huán)境下溫度變化較小， 通常將材料熱工參數(shù)視為常數(shù)， 考慮了3種熱傳遞方式對(duì)結(jié)構(gòu)溫度場(chǎng)的影響； 李國(guó)強(qiáng)等[18]在對(duì)火災(zāi)下二維鋼板-混凝土組合樓板溫度場(chǎng)分析時(shí)， 考慮結(jié)構(gòu)的熱對(duì)流和熱輻射邊界條件， 并探究火災(zāi)位置和分布對(duì)樓板溫度場(chǎng)的影響， 通過(guò)對(duì)結(jié)果進(jìn)行參數(shù)分析， 提出了適用于火災(zāi)下的樓板溫度計(jì)算的簡(jiǎn)化公式； 胡克旭等[19]采用有限單元法， 假定單元溫度隨著位置坐標(biāo)線性變化， 計(jì)算了火災(zāi)下壓型鋼板-混凝土組合樓板在存在熱對(duì)流和熱輻射， 同時(shí)考慮建筑材料熱工參數(shù)隨溫度變化時(shí)的溫度場(chǎng).

實(shí)際上， 由于材料比熱和熱傳導(dǎo)率均會(huì)隨著溫度變化， 加之邊界條件的復(fù)雜性， 使得求解偏微分方程得出解析解是極少數(shù)情況. 采用數(shù)值方法計(jì)算， 能有效簡(jiǎn)化計(jì)算過(guò)程. 常見的數(shù)值解法有有限元法、有限體積法、邊界元法和有限差分法等， 這些方法中， 有限差分法具有簡(jiǎn)單且穩(wěn)定的優(yōu)點(diǎn)[20]. 本研究基于熱傳導(dǎo)的有限差分方程， 使用MATLAB編程計(jì)算結(jié)構(gòu)溫度場(chǎng). 以簡(jiǎn)化的一維變系數(shù)瞬態(tài)溫度場(chǎng)分布模型為例， 闡明考慮材料熱工參數(shù)隨溫度變化時(shí)的有限差分法計(jì)算過(guò)程， 并建立有限元模型進(jìn)行結(jié)果驗(yàn)證.

1 熱傳導(dǎo)差分方程的建立

理論推導(dǎo)求解物體溫度場(chǎng)一般采用解析法和數(shù)值法. 其中， 解析法是以熱傳導(dǎo)方程為基礎(chǔ)， 求解偏微分方程， 尋求物體溫度分布與時(shí)間和位置關(guān)系的函數(shù)表達(dá)式， 即T(x，y，z，t). 實(shí)際能采用解析法求解的熱傳導(dǎo)問(wèn)題是極少數(shù)情況， 因此， 通常采用有限差分的數(shù)值方法建立傳熱的有限差分方程， 求解溫度場(chǎng)近似值[21]， 通過(guò)減小單元格大小(時(shí)間步長(zhǎng)、空間步長(zhǎng))來(lái)提高計(jì)算精度.

有限差分方程的建立可通過(guò)兩種方式完成： 一是基于熱傳導(dǎo)方程建立差分方程， 二是基于能量平衡建立差分方程.

1.1 基于熱傳導(dǎo)方程的差分方程

熱傳導(dǎo)方程為：

(1)

式中：λ(T)為熱傳導(dǎo)率；ρ為材料密度；C(T)為比熱容.(1)式表達(dá)了溫度在空間上的分布與溫度隨時(shí)間變化的關(guān)系：

基于熱傳導(dǎo)方程建立傳熱的差分方程， 就是以熱傳導(dǎo)方程為基礎(chǔ)， 用差分代替微分， 將偏導(dǎo)寫成差分形式進(jìn)行求解. 對(duì)于任意連續(xù)函數(shù)， 使用泰勒級(jí)數(shù)將x0相鄰兩點(diǎn)函數(shù)值f(x-1)和f(x1)在x0處展開得：

(2)

(3)

(2)式+(3)式得到f(x)在x0處沿著x方向二階導(dǎo)數(shù)， 略去高階項(xiàng)， 得到：

(4)

同理， 在圖1所示三維直角坐標(biāo)系中， 與點(diǎn)P(x0，y0，z0)相鄰的上下左右前后6個(gè)點(diǎn)的溫度值分別為TU，TD，TL，TR，TF，TB.

圖1 空間位置示意

則T(x，y，z)在點(diǎn)P(x0，y0，z0)處各方向二階偏導(dǎo)可表示為：

(5)

取空間間隔Δx=Δy=Δz=Δ， 時(shí)間間隔為Δt， 時(shí)間偏導(dǎo)采用向前差分格式則有：

(6)

式中：T′為P點(diǎn)下一時(shí)刻溫度值.將(5)式和(6)式帶入(1)式中， 得到基于熱傳導(dǎo)方程的差分方程， 整理得：

(7)

(8)

1.2 基于能量平衡的差分方程

基于能量平衡建立差分方程， 是根據(jù)單元能量平衡的思想進(jìn)行的， 即： 單元內(nèi)輸入熱量和輸出熱量的差值， 引起單元溫度變化. 差值為正， 單元溫度升高； 差值為負(fù)， 單元溫度降低.

在圖2所示三維空間直角坐標(biāo)系中取控制單元(圖3).

圖2 空間直角坐標(biāo)

圖3 控制單元示意

圖2中A為外界節(jié)點(diǎn)(邊界條件)，P1為邊界節(jié)點(diǎn)，P2為中間節(jié)點(diǎn). 圖3中， 與點(diǎn)P相鄰的6個(gè)點(diǎn)U，D，L，R，F(xiàn)，B分別代表點(diǎn)P上下左右前后相鄰點(diǎn). 為簡(jiǎn)化計(jì)算， 令網(wǎng)格間隔Δx=Δy=Δz=Δ， 因此相鄰各點(diǎn)到P的距離均為Δ. 節(jié)點(diǎn)i到相鄰節(jié)點(diǎn)j之間的傳熱系數(shù)為Kij， 基于傅里葉定律， 從節(jié)點(diǎn)j到i的熱量為[13]：

Qij=Kij(Tj-Ti)

(9)

單元內(nèi)蓄積的能量使得單元內(nèi)部溫度發(fā)生變化， 節(jié)點(diǎn)i的能量平衡方程為：

(10)

圖3中， 中間節(jié)點(diǎn)

(11)

圖2中， 邊界節(jié)點(diǎn)(假設(shè)為左側(cè)邊界， 右側(cè)同理)

(12)

(13)

對(duì)(12)式右側(cè)進(jìn)行向前差分， 左側(cè)采用加權(quán)平均的形式， 整理得：

(14)

式中， 若β=0， 方程為顯式格式；β=1， 方程為隱式格式. 得到差分方程， 結(jié)合相應(yīng)的邊界條件即可得到溫度場(chǎng)分布.

1.3 差分方程的特殊形式

實(shí)際情況中多為三維熱傳導(dǎo)問(wèn)題. 但在處理一些特定問(wèn)題時(shí)， 適當(dāng)做出假設(shè)， 可將三維簡(jiǎn)化為二維或一維情況.

按照1.2節(jié)思路， 可得二維熱傳導(dǎo)差分方程為：

(15)

一維熱傳導(dǎo)差分方程為：

(16)

同于1.1節(jié)， 若β=0， 方程為顯式格式；β=1， 方程為隱式格式.

1.4 顯示差分的穩(wěn)定性要求

對(duì)于顯式格式的有限差分方程， 在計(jì)算時(shí)， 涉及到方程的穩(wěn)定性問(wèn)題.

令式(14)中β=0， 得顯式格式差分方程：

(17)

由(17)式可知， 由于FO>0， 所以當(dāng)T的系數(shù)(1/FO-6)<0時(shí)，T越大T′ 就會(huì)越小， 即此時(shí)刻溫度越大下一時(shí)刻的溫度越小， 顯然違背熱力學(xué)原理. 因此， (1/FO-6)≥0同理. 對(duì)于二維和一維有相同的約束.

內(nèi)部節(jié)點(diǎn)：

一維直角坐標(biāo)FO≤1/2

二維直角坐標(biāo)FO≤1/4

三維直角坐標(biāo)FO≤1/6

(18)

邊界節(jié)點(diǎn)(對(duì)流邊界條件)：

一維直角坐標(biāo)FO≤1/[2(1+Bi)]

二維直角坐標(biāo)FO≤1/[2(2+Bi)]

三維直角坐標(biāo)FO≤1/[2(3+Bi)]

(19)

式中：Bi=hΔ/λ， 時(shí)間步長(zhǎng)限制取(18)式與(19)式所得交集. 在一維直角坐標(biāo)系中， 由(19)式可得時(shí)間步長(zhǎng)限制為：

(20)

2 差分方程求解

以“一維無(wú)熱源具有熱對(duì)流邊界條件的熱傳導(dǎo)問(wèn)題”為例， 分別介紹是否考慮材料熱工參數(shù)隨溫度變化兩種情況時(shí)的溫度場(chǎng)求解方法.

已知初始時(shí)刻所有節(jié)點(diǎn)溫度均為T(X， 0)=T0； 左邊界環(huán)境溫度滿足TA=A(t)， 右邊界環(huán)境溫度滿足TB=B(t)， 表面?zhèn)鳠嵯禂?shù)為hw. 將構(gòu)件沿x方向均勻劃分為N個(gè)單元格(圖4).

圖4 單元?jiǎng)澐质疽?/p>

2.1 常系數(shù)熱傳導(dǎo)方程

在上述一維熱傳導(dǎo)問(wèn)題中， 取T(X， 0)=T0=20 ℃， 左側(cè)環(huán)境溫度取國(guó)際標(biāo)準(zhǔn)升溫曲線(ISO834)[14]溫度值A(chǔ)(t)=20+345lg(8t+1)，B(t)=20 ℃，L=4 m，N=80. 熱工參數(shù)見表1.

表1 熱傳導(dǎo)相關(guān)參數(shù)值

對(duì)于常系數(shù)熱傳導(dǎo)問(wèn)題， 采用顯式和隱式格式均可求解. 由(13)式和(16)式得節(jié)點(diǎn)i的差分方程， 其中顯式差分格式為：

(21)

隱式差分格式為：

(22)

2.1.1 常系數(shù)顯式求解

采用顯式格式求解熱傳導(dǎo)問(wèn)題時(shí)， 首先根據(jù)精度要求， 選取空間網(wǎng)格Δ大小， 再根據(jù)顯式格式的穩(wěn)定性要求確定時(shí)間步長(zhǎng)Δt， 最后依據(jù)遞推關(guān)系(21)式， 結(jié)合邊界和初始條件計(jì)算得到溫度場(chǎng).

Δt≤49.062 5 s

計(jì)算時(shí)間取100 min， 則所需時(shí)間步為：

m=122.29

取Δt=49.062 5 s時(shí)，F(xiàn)O=0.187 5，Bi=5/3， 帶入(21)式， 得：

(23)

至此得到該例的顯式差分格式， 已知初始溫度分布T(x， 0)和邊界條件TA，TB， 依次沿著時(shí)間步向下計(jì)算， 可得到所有節(jié)點(diǎn)任意時(shí)刻溫度值. 計(jì)算示意圖見圖5， 橫軸為節(jié)點(diǎn)位置，i表示第i個(gè)節(jié)點(diǎn)； 縱軸為時(shí)間步，j表示第j個(gè)時(shí)間步.

圖5 顯式差分示意圖

由圖5可見， 顯示格式遞推關(guān)系清晰， 易于計(jì)算， 根據(jù)(23)式采用MATLAB編程計(jì)算溫度場(chǎng)即可.

2.1.2 常系數(shù)隱式求解

采用隱式格式求解熱傳導(dǎo)問(wèn)題時(shí)， 根據(jù)精度要求， 選取空間網(wǎng)格Δ大小， 隱式格式的差分方程沒(méi)有穩(wěn)定性要求， 時(shí)間間隔可根據(jù)需求選定. 通過(guò)(22)式， 結(jié)合邊界和初始條件， 求解方程組， 得到溫度場(chǎng).

初始時(shí)刻， 時(shí)間步j(luò)=0. 由(22)式得每個(gè)節(jié)點(diǎn)上的差分方程為：

寫成矩陣形式：

(24)

其中：K1=1/(2FO)+(1+Bi)·2FO；K=2FO·Bi.

同理， 對(duì)于下一個(gè)時(shí)間步， 也可以得到相似的矩陣形式. 由此可見： 對(duì)于隱式格式， 雖然方程的計(jì)算不受時(shí)間間隔的影響， 但是在求解每一時(shí)刻各節(jié)點(diǎn)溫度值時(shí)， 都需要聯(lián)立方程組， 這使得求解計(jì)算量變大. 隱式格式求解結(jié)構(gòu)圖見圖6.

圖6 隱式差分示意圖

此時(shí)基于矩陣(24)求解， 所采用的方法稱為三對(duì)角矩陣算法(TDMA). 取時(shí)間間隔為20 s， 空間間隔為0.05 m， 采用MATLAB編程求解矩陣， 即可得到溫度場(chǎng).

2.2 變系數(shù)熱傳導(dǎo)方程

實(shí)際上， 當(dāng)溫度較高時(shí)， 材料熱工參數(shù)會(huì)隨著溫度的改變而改變， 即熱傳導(dǎo)系數(shù)和比熱為溫度的函數(shù). 在計(jì)算考慮材料熱工參數(shù)隨著溫度的溫度場(chǎng)時(shí)， 本研究采用歐洲規(guī)范EC3[2]中鋼材熱傳導(dǎo)系數(shù)和比熱容隨溫度變化關(guān)系， 由于材料熱工參數(shù)為溫度的函數(shù)， 熱傳導(dǎo)方程變?yōu)榉蔷€性方程， 基于熱傳導(dǎo)方程直接求解會(huì)十分復(fù)雜. 采用能量平衡的方法建立顯式格式的差分方程， 計(jì)算時(shí)遵循以下流程： 在已知初始溫度場(chǎng)T0時(shí)， 依據(jù)EC3規(guī)范計(jì)算得到初始時(shí)刻的熱工參數(shù)λ0和C0， 再通過(guò)λ0和C0依據(jù)差分方程計(jì)算新的溫度場(chǎng)T1. 即已知前一時(shí)刻熱工參數(shù)， 可得到下一時(shí)刻新的溫度場(chǎng)， 多次計(jì)算得任意時(shí)刻的溫度場(chǎng)Tn. 依據(jù)以上思路和1.2節(jié)內(nèi)容， 建立差分方程.

對(duì)于中間節(jié)點(diǎn)， 有：

(25)

(26)

對(duì)于左側(cè)邊界節(jié)點(diǎn)， 有：

(27)

取顯式格式：

(28)

同理， 右側(cè)邊界節(jié)點(diǎn)：

(29)

若仍取Δ=0.05 m. 由于采用顯式差分格式， 需滿足(20)式與(22)式的穩(wěn)定性要求. 結(jié)合(25)式和(26)式可知， 取初始的熱工參數(shù)λ0和C0滿足要求， 帶入?yún)?shù)得：

Δt≤21.36 s 取Δt=20 s

可以看出， 方程中溫度的系數(shù)都是隨溫度變化的， 即非線性方程， 采用2.1節(jié)中直接求解的方法無(wú)法得到結(jié)果. 根據(jù)(28)式-(29)式， 考慮溫度引起的材料性能改變(遵循EC3規(guī)范公式)， 使用MATLAB編程即可計(jì)算得到溫度場(chǎng).

3 有限元對(duì)比與分析

基于以上內(nèi)容， 使用ABAQUS分別建立不考慮熱工參數(shù)隨溫度變化的有限元模型MODEL0與考慮熱工參數(shù)隨溫度變化的有限元模型MODEL1， 相關(guān)參數(shù)與2.1節(jié)、2.2節(jié)保持一致. 模型MODEL0與模型MODEL1均采用C3D8T單元， 分別取x=0 m和x=1 m處截面溫度進(jìn)行對(duì)比.

3.1 常系數(shù)計(jì)算結(jié)果驗(yàn)證

不考慮材料性能隨溫度變化時(shí)， 可采用顯式和隱式兩種差分格式計(jì)算， 將2.1節(jié)計(jì)算結(jié)果與MODEL0計(jì)算結(jié)果進(jìn)行對(duì)比， 以驗(yàn)證計(jì)算結(jié)果的準(zhǔn)確性(圖7和圖8).

由圖7可知， 當(dāng)不考慮材料熱工參數(shù)隨溫度變化時(shí)， 采用隱式差分計(jì)算結(jié)果與ABAQUS有限元模擬結(jié)果高度吻合. 在環(huán)境溫度采用國(guó)際標(biāo)準(zhǔn)升溫曲線(ISO834)時(shí)， 在邊界處， 物體內(nèi)部溫度變化與環(huán)境溫度變化形式相同， 均呈對(duì)數(shù)形式； 在遠(yuǎn)離邊界的物體內(nèi)部， 溫度變化幅度隨著離邊界距離的增大而減小， 且溫度變化呈指數(shù)形式. 對(duì)比圖8可知， 采用顯式差分計(jì)算在x=0 m邊界處與ABAQUS有限元模擬結(jié)果幾乎重合； 在距離邊界1 m處， 顯式計(jì)算結(jié)果和ABAQUS模擬結(jié)果隨著時(shí)間的增加相差變大， 但總體差距仍然較小(小于0.03 ℃). 產(chǎn)生這一現(xiàn)象的原因在于， 隱式差分計(jì)算中采用的TDMA方法計(jì)算結(jié)果是理論上的精確解， 而顯式差分的精度取決于網(wǎng)格劃分的大小.

圖7 常系數(shù)隱式計(jì)算結(jié)果

圖8 常系數(shù)顯式計(jì)算結(jié)果

3.2 變系數(shù)計(jì)算結(jié)果驗(yàn)證

當(dāng)考慮材料性能隨溫度變化時(shí)， 采用顯式差分格式計(jì)算， 將2.2節(jié)計(jì)算結(jié)果與MODEL1進(jìn)行對(duì)比， 以驗(yàn)證計(jì)算結(jié)果的準(zhǔn)確性(圖9)所示.

圖9 變系數(shù)計(jì)算結(jié)果

由圖9可見， 當(dāng)考慮材料熱工參數(shù)隨溫度變化時(shí)， 使用2.2節(jié)中方法計(jì)算結(jié)果與ABAQUS有限元模擬結(jié)果吻合度較高. 在環(huán)境溫度采用國(guó)際標(biāo)準(zhǔn)升溫曲線， 考慮材料熱工參數(shù)隨溫度變化時(shí)， 在邊界處物體內(nèi)部溫度變化與環(huán)境溫度變化形式相同， 呈對(duì)數(shù)形式； 遠(yuǎn)離邊界的物體內(nèi)部， 溫度變化范圍隨著離邊界距離的增大而減小， 溫度變化呈指數(shù)形式. 溫度場(chǎng)變化規(guī)律與常系數(shù)相同. 可以看到， 采用有限差分法計(jì)算時(shí)， 不論是常系數(shù)還是變系數(shù)計(jì)算結(jié)果與相應(yīng)的有限元模擬結(jié)果均存在誤差， 這些誤差主要來(lái)源于兩方面： ①采用有限元軟件計(jì)算的溫度場(chǎng)實(shí)際上是采用有限單元法計(jì)算的， 本研究采用的是有限差分法， 不論是有限單元法還是有限差分法， 兩者作為數(shù)值方法本身就與精確解存在誤差； ②兩種方法的計(jì)算精度不僅與時(shí)間網(wǎng)格、空間網(wǎng)格大小的有關(guān)， 有限單元法中選取的形函數(shù)、有限差分法中省略的高階項(xiàng)也直接影響相應(yīng)的計(jì)算精度.

3.3 常系數(shù)、變系數(shù)計(jì)算結(jié)果對(duì)比

為探究是否考慮材料熱工參數(shù)對(duì)溫度場(chǎng)的影響， 分別取離左側(cè)邊界x=0 m，x=1 m處截面兩種情況下顯示格式計(jì)算所得溫度場(chǎng)進(jìn)行對(duì)比， 并進(jìn)行誤差分析， 結(jié)果見圖10和圖11.

圖10 變系數(shù)和常系數(shù)溫度場(chǎng)對(duì)比

圖11 變系數(shù)和常系數(shù)溫度場(chǎng)誤差對(duì)比

由圖10可見， 是否考慮材料熱工參數(shù)隨溫度變化不影響溫度分布規(guī)律， 但改變了溫度大小： 考慮材料熱工參數(shù)隨溫度變化時(shí)的溫度相對(duì)于不考慮材料熱工參數(shù)隨溫度變化時(shí)的溫度偏高. 由圖11可見： 采用有限差分法時(shí)不論是否考慮材料熱工參數(shù)隨溫度變化， 兩者所產(chǎn)生的誤差具有相同的變化規(guī)律， 但不考慮材料熱工參數(shù)隨溫度變化產(chǎn)生的誤差總是高于考慮熱工參數(shù)隨溫度變化的誤差； 圖11a中， 除時(shí)間t=20 s時(shí)誤差為60%， 其余各個(gè)時(shí)間誤差均小于10%， 這是由于采用的升溫曲線初始溫度變化快， 同時(shí)計(jì)算時(shí)所取時(shí)間步長(zhǎng)較大引起的. 因此， 當(dāng)環(huán)境溫度較高或溫度變化速率較大時(shí)， 計(jì)算物體溫度場(chǎng)時(shí)應(yīng)考慮材料熱工參數(shù)隨溫度的變化， 且取更小的時(shí)間空間網(wǎng)格以提高精度.

4 結(jié)論

本研究在推導(dǎo)求解變系數(shù)熱傳導(dǎo)問(wèn)題的基礎(chǔ)上， 以具有熱對(duì)流邊界條件同時(shí)考慮材料熱工參數(shù)隨溫度變化的一維瞬態(tài)熱傳導(dǎo)問(wèn)題為例， 采用有限差分法分析了當(dāng)環(huán)境溫度為國(guó)際標(biāo)準(zhǔn)升溫曲線時(shí)物體內(nèi)部溫度分布情況， 得到以下五點(diǎn)結(jié)論.

1) 基于能量平衡建立傳熱有限差分方程相對(duì)于從熱傳導(dǎo)方程出發(fā)建立傳熱有限差分方程， 過(guò)程清晰、可考慮熱工參數(shù)隨溫度變化的情況， 更適用于特殊復(fù)雜情況.

2) 通過(guò)與相應(yīng)的有限元模型計(jì)算結(jié)果對(duì)比表明： 本研究采用基于能量平衡建立的有限差分解法能有效計(jì)算考慮材料熱工參數(shù)隨溫度變化時(shí)的溫度場(chǎng).

3) 是否考慮材料熱工參數(shù)隨溫度變化不影響溫度分布規(guī)律， 但影響溫度大小： 考慮材料熱工參數(shù)隨溫度變化時(shí)得到的溫度場(chǎng)， 高于不考慮熱工參數(shù)變化得到的溫度場(chǎng).

4) 在環(huán)境溫度采用國(guó)際標(biāo)準(zhǔn)升溫曲線(ISO834)時(shí)， 在邊界處， 物體內(nèi)部溫度變化與環(huán)境溫度變化形式相同， 均呈對(duì)數(shù)形式； 在遠(yuǎn)離邊界的物體內(nèi)部， 溫度變化幅度隨著離邊界的距離增大而減小， 且溫度變化呈指數(shù)形式.

5) 依據(jù)本研究的思路也可計(jì)算二維、三維溫度場(chǎng)， 進(jìn)一步就特定的環(huán)境溫度提出相應(yīng)的物體內(nèi)部溫度場(chǎng)簡(jiǎn)化計(jì)算方法.