側面作用移動熱源柱體溫度場適體坐標有限容積法數值解與分析解的對比

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

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

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

1 問題的物理數學描述

1.1 有量綱參數描述

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

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

問題的數學描述為能量守恒定律：

(1)

方程(1)是非穩態熱傳導微分方程,其中ρ為密度,cp為比熱容,T為溫度,t為時間,為拉普拉斯算子,λ為導熱系數.

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

熱源作用區2φ0作用的熱通量

(2)

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

(3)

這里h是均勻的.

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

T0=Tf=0,t≤0.

(4)

1.2 無量綱參數描述

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

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

(5)

(6)

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

(7)

則上述描述變為：

(8)

在熱源作用區

(9)

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

(10)

初始條件

θ=θ0,t*≤0.

(11)

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

(12)

2 數學模型的解析解

對于上述問題文獻[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的數量級,從而即便是中等角速度,ωn也有較大的數值,根據I(x)的特性,I(ωn)在較小的n條件下就會有I(ωn)→∞.這樣會使得式(15)～(16)的計算沒有意義(數值為∞).在這種情況下需要把bern(x)及bein(x)表示成另外一種形式,使得M1(x),M2(x),M3(x)成c1·exp(c2x)的形式,然后帶入式(14),消去exp(x)后就可保證式(14)的數值有意義(不會是∞).

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

(26)

(27)

再由遞推式(21)有：

(28)

(29)

(30)

(31)

(32)

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

(33)

(34)

(35)

3 數學模型的數值解方法

3.1 數值計算區域

用數值方法求解以上問題時選取整體區域并劃分網格為如圖2所示.

圖2 整體計算時的網格Fig.2 Grid of total calculation

3.2 方程的變換[16]

將物理空間用適體坐標系轉化到計算空間,如圖3所示.

圖3 適體坐標系的轉換及控制容積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)代數方程組及6個面上的邊界代數方程.求解這一方程組,就可獲得溫度場.論文采用TDMA方法求解代數方程組.收斂判據是同一時層的迭代過程中最大溫度相對誤差不超過10-5.

圖4 控制容積Fig.4 control volume

3.4 移動熱源的處理方法

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

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

4 結果及分析

4.1 計算用模型說明

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

4.2 解析解的特性

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

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

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

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

表1 特征值μ0,j的值

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

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

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

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

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

4.3 數值解與解析解的比較

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

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

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

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

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

圖11 同一時刻圓柱側表面上無量綱溫度隨位置變化的解析解與數值解對比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 結論

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

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

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

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