側(cè)面作用移動熱源柱體溫度場適體坐標(biāo)有限容積法數(shù)值解與分析解的對比

王 璐，石廣田，翟婉明,3，王良璧

(1. 蘭州交通大學(xué) 機電工程學(xué)院，蘭州 730070；2. 鐵道車輛熱工教育部重點實驗室(蘭州交通大學(xué))，蘭州 730070；3. 西南交通大學(xué)牽引動力國家重點實驗室，成都 610031)

由摩擦生熱產(chǎn)生的移動熱源會對車輪及鋼軌的性能產(chǎn)生較大的影響[1-7].從而,對車輪和鋼軌在移動熱源的作用下溫度場特性的研究非常重要[8-10].由于這類移動熱源問題的復(fù)雜性無法用解析法獲得溫度場.又由于鐵道車輛輪軌工作的特殊性,一比一的實驗方法也非常難于實施,實驗室中的實驗數(shù)據(jù)由于很難轉(zhuǎn)化到可應(yīng)用到實際情況的數(shù)據(jù)而受到了限制.這樣數(shù)值方法成為研究車輪在移動熱源作用下溫度場特性的有效方法.從而非常有必要發(fā)展該類問題的數(shù)值方法及實施的程序.目前在用的商業(yè)軟件很多,但大多數(shù)軟件處理有追趕特點移動熱源問題的邊界條件時非常復(fù)雜[10].另外,關(guān)于數(shù)值方法的嚴(yán)格考核報道的較少.鑒于這些原因,論文擬發(fā)展一適合處理追趕特點移動熱源邊界的有限容積法數(shù)值分析程序.并與分析解進行比對,驗證其數(shù)值分析程序的正確性.

1 問題的物理數(shù)學(xué)描述

1.1 有量綱參數(shù)描述

圖1所示是半徑為r0的一無限長柱體,在柱體側(cè)面上作用有一段AB的移動熱源,移動速度為v,移動熱源的強度為均勻強度qb.作用面積是2φ0對應(yīng)的弧長和柱體長度的乘積.除移動熱源加熱外,在柱體表面上有表面?zhèn)鳠嵯禂?shù)為常數(shù)h的對流傳熱過程.物理問題是用數(shù)值方法確定其溫度場.

圖1 側(cè)面作用移動熱源的無限長柱體模型Fig.1 Infinite cylinder model of moving heat source action on the side surface

問題的數(shù)學(xué)描述為能量守恒定律：

(1)

方程(1)是非穩(wěn)態(tài)熱傳導(dǎo)微分方程,其中ρ為密度,cp為比熱容,T為溫度,t為時間,為拉普拉斯算子,λ為導(dǎo)熱系數(shù).

方程(1)的邊界條件為：

熱源作用區(qū)2φ0作用的熱通量

(2)

側(cè)面對流傳熱區(qū)(包含熱源作用區(qū)與非作用區(qū))

(3)

這里h是均勻的.

方程(1)的初始條件為：

T0=Tf=0,t≤0.

(4)

1.2 無量綱參數(shù)描述

為了無量綱化方程(1),這里無量綱定義如下：

θ=πλvT/(2αqb).

(5)

(6)

X=vx/2α,Y=vy/2α,Z=vz/2α.

(7)

則上述描述變?yōu)椋?/p>

(8)

在熱源作用區(qū)

(9)

側(cè)面對流傳熱區(qū)(包含熱源作用區(qū)與非作用區(qū))

(10)

初始條件

θ=θ0,t*≤0.

(11)

θw=πλvTw/(2αqb),θf=0,θ0=0.

(12)

2 數(shù)學(xué)模型的解析解

對于上述問題文獻[14]中給出了解析解：

(13)

這里

(14)

及

(15)

(16)

這里μn,j是下列式子的根.

(17)

這里有：

(18)

J-n(x)=(-1)nJn(x).

(19)

(20)

(21)

(22)

bernx+ibeinx=inIn(xi1/2).

(23)

(24)

(25)

對于鋼件,由于α比較小,大約是10-6的數(shù)量級,從而即便是中等角速度,ωn也有較大的數(shù)值,根據(jù)I(x)的特性,I(ωn)在較小的n條件下就會有I(ωn)→∞.這樣會使得式(15)～(16)的計算沒有意義(數(shù)值為∞).在這種情況下需要把bern(x)及bein(x)表示成另外一種形式,使得M1(x),M2(x),M3(x)成c1·exp(c2x)的形式,然后帶入式(14),消去exp(x)后就可保證式(14)的數(shù)值有意義(不會是∞).

文獻[15]給出了bern(x)及bein(x)這種形式：

(26)

(27)

再由遞推式(21)有：

(28)

(29)

(30)

(31)

(32)

把式(26)～(32)代入式(14)就有θperiodic近似式：

(33)

(34)

(35)

3 數(shù)學(xué)模型的數(shù)值解方法

3.1 數(shù)值計算區(qū)域

用數(shù)值方法求解以上問題時選取整體區(qū)域并劃分網(wǎng)格為如圖2所示.

圖2 整體計算時的網(wǎng)格Fig.2 Grid of total calculation

3.2 方程的變換[16]

將物理空間用適體坐標(biāo)系轉(zhuǎn)化到計算空間,如圖3所示.

圖3 適體坐標(biāo)系的轉(zhuǎn)換及控制容積Fig.3 Transformation and control volume of Body-Fitted coordinate system

令：

(36)

在計算空間方程(8)變換為：

(37)

3.3 方程的離散

以上方程可以圖4所示的控制容積中離散.

aPθP=aEθE+aWθW+aNθN+aSθS+aTθT+aBθB+b.

(38)

方程離散后獲得一個(L1-2)×(M1-2)×(N1-2)代數(shù)方程組及6個面上的邊界代數(shù)方程.求解這一方程組,就可獲得溫度場.論文采用TDMA方法求解代數(shù)方程組.收斂判據(jù)是同一時層的迭代過程中最大溫度相對誤差不超過10-5.

圖4 控制容積Fig.4 control volume

3.4 移動熱源的處理方法

移動熱源在數(shù)值計算過程中需特殊處理.如圖5所示,把柱體側(cè)面展開,AB線是同一條線,這時熱源分布在AB周圍及熱源分布在兩AB線之內(nèi).為保證計算過程中加入柱體的熱量守恒,沿移動方向,每一時間步長熱源移動一完整網(wǎng)格.

圖5 計算區(qū)域移動熱源分布示意圖Fig.5 Distribution diagram of moving heat source in calculation area

4 結(jié)果及分析

4.1 計算用模型說明

圖1模型和圖2網(wǎng)格所示的圓柱半徑為a=0.152 4 m.熱源作用區(qū)為2φ0=π/10 rad,計算使用時間步長為Δt=0.001 s,α=4.415 955×10-6m2/s,λ=15.24 W/m·K,ω=10π,轉(zhuǎn)一圈時間周期=2π/ω=0.2 s,無量綱時間周期為0.2α/a2=0.2×4.415 955×10-6/0.152 42=3.802 635×10-6.一個點從熱源進入到退出所持續(xù)時間為φ0/ω=1/100 s.無量綱時間步長Δt*=0.001×4.416×10-6/0.152 42=1.901×10-7.圓柱表面對流傳熱系數(shù)h=1 000 W/m2·K,畢沃?jǐn)?shù)NB=10.

4.2 解析解的特性

方程(13)第一項的特性如圖6所示,其值隨時間的變化規(guī)律為周期性變化.但其幅值不隨時間變化.為了表明這一特性,我們把解析解第一項很長時間內(nèi)θperiodic(周期性值)的變化如圖6(a)所示.為了減小該圖的字節(jié)數(shù),圖中只畫出了三個較小時間段的溫度波線.可以看出,θperiodic在較長時間段內(nèi)的峰值變化并不明顯.每一段的波線在圖6(b)中放大.有圖6(b)可知,波的幅值隨時間變化非常小.

圖6 解析解第一項(周期性值)隨時間的變化規(guī)律Fig.6 The changes with time of first term (periodic value) of the analytic solution

方程(13)第二項(用θ2表示)的特性如圖7所示,其數(shù)值和j的大小有關(guān).方程(17)n=0的特征值μ0,j如表1所示.總體而言,無論j取何值,θ2的變化規(guī)律都是隨時間的增加而逐漸增大的.隨j的增大,θ2的初始值越小,增大的速度越快.當(dāng)j>50時無量綱溫度約接近數(shù)值解.

圖7 解析解第二項(周期性值)隨時間的變化規(guī)律Fig.7 The changes with time of second term (periodic value) of the analytic solution

表1 特征值μ0,j的值

方程(13)第三項(用θ3表示)的特性如圖8所示,其數(shù)值和j的大小有關(guān).但相比θperiodic和θ2兩項,其數(shù)量級較小,且僅在移動熱源剛開始作用時刻有微小作用,隨時間的增加,第三項數(shù)值減小為接近0,對方程(13)的幾乎沒有影響.

圖8 解析解第三項隨時間的變化規(guī)律Fig.8 The changes with time of third term of the analytic solution

解析解獲得的圓柱側(cè)表面一點無量綱溫度隨時間變化規(guī)律如圖9所示.

圖9 圓柱側(cè)表面一點無量綱溫度隨時間的變化規(guī)律Fig.9 Dimensionless temperature changing with the time of one point on the cylindrical side surface

無量綱溫度在熱源剛剛作用時躍升,熱源剛離開的一段時間內(nèi)驟降,溫度下降到在較接近環(huán)境溫度后,下降趨勢減緩,熱源再次作用時,無量綱溫度再次上升.由于周期性運動,無量綱溫度也呈現(xiàn)周期性規(guī)律.在無限長圓柱剛開始轉(zhuǎn)動的時候,各個周期的最高溫度隨時間有升高的趨勢,在運動一段時間后,熱源造成的最大溫升趨于穩(wěn)定,可以看到側(cè)表面溫度隨時間的變化為基本穩(wěn)定的周期性規(guī)律,最高溫度和最低溫度基本為定值.

4.3 數(shù)值解與解析解的比較

1) 圓柱側(cè)表面一點無量綱溫度隨時間變化規(guī)律的解析解與數(shù)值解對比

圓柱側(cè)表面一點無量綱溫度隨時間變化規(guī)律的解析解與數(shù)值解對比如圖10所示,無論是采用解析解還是數(shù)值解,得到的無量綱溫度都具有相同的隨時間變化的規(guī)律,且對比各個時刻的無量綱溫度,兩種方法的計算結(jié)果誤差非常小.可以驗證數(shù)值方法的準(zhǔn)確性.

2) 同一時刻圓柱側(cè)面無量綱溫度分布的解析解與數(shù)值解比較

同一時刻圓柱側(cè)面無量綱溫度分布隨位置變化的解析解與數(shù)值解比較如圖11所示,可以看到,在同一時刻不同位置的無量綱溫度存在著相似的變化規(guī)律.在熱源作用區(qū),無量綱溫度突然上升并在即將離開熱源作用區(qū)的位置達到最大值,越遠離熱源作用區(qū),無量綱溫度越低,其值減小的趨勢越來平緩.但無論是采用解析解還是數(shù)值解,得到的無量綱溫度都具有相同的隨位置變化的規(guī)律,且對比同一時刻各個位置的無量綱溫度,兩種方法的計算結(jié)果誤差非常小,驗證了數(shù)值方法的準(zhǔn)確性.

圖10 圓柱側(cè)表面一點溫度隨時間變化的解析解與數(shù)值解對比Fig.10 A comparison between the analytical results and numerical results of the temperature variation with time on the cylinder side surface

圖11 同一時刻圓柱側(cè)表面上無量綱溫度隨位置變化的解析解與數(shù)值解對比Fig.11 A comparison between the analytical results and numerical results of the dimensionless temperature variation on the cylinder side surface varies with position at the same time

5 結(jié)論

為了獲得鐵道車輛車輪在接近實際移動熱源作用下的溫度場,論文發(fā)展了一有限容積法數(shù)值方法程序,對該程序的準(zhǔn)確性進行了驗證.驗證方法是用數(shù)值方法和解析方法解一無限長圓柱受一移動熱源作用的溫度場,比較兩結(jié)果的差異.獲得的結(jié)論如下:

1) 論文發(fā)展的采用有限容積法的程序獲得的數(shù)值解和對應(yīng)分析解相差甚小,程序計算結(jié)果準(zhǔn)確性高,這說明發(fā)展的程序可以用來進行移動熱源引起溫度場的計算.

2) 解析解中第一項是隨時間周期變化的,但其變化的幅值幾乎不隨時間的變化而變化.解析中第二項隨時間是漸進增加的.第三項較第一及第二項小,在移動熱源開始作用時刻有微小作用,隨時間的增加,第三項數(shù)值減小為接近0.

3) 獲得分析解的計算工作量不亞于數(shù)值方法的計算工作量,這與常規(guī)的認(rèn)識有較大的出入.