熱源周期振蕩條件下一維融化問題的數(shù)值研究

曲良輝，楊金根，邢 琳，令 鋒，董 麗

(1. 中原工學(xué)院 理學(xué)院，河南 鄭州 450007；2. 信陽(yáng)師范學(xué)院 數(shù)學(xué)與信息科學(xué)學(xué)院，河南 信陽(yáng) 464000；3. 肇慶學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，廣東 肇慶 526061)

0 引言

在自然界和各種工業(yè)生產(chǎn)中，經(jīng)常會(huì)遇到這些過程，如冰的融化、金屬的重結(jié)晶、合金的凝固、食物的冷藏、液滴的蒸發(fā)和粒子的擴(kuò)散等，它們都涉及相變傳熱問題[1-3]．其特點(diǎn)是區(qū)域內(nèi)存在一個(gè)隨時(shí)間變化的兩相界面，在該界面上吸收或放出熱量．這類問題通常被稱為Stefan問題，也稱為移動(dòng)邊界問題，因?yàn)閮上嘟缑娴奈恢檬且苿?dòng)的．在數(shù)學(xué)上相變傳熱問題是一個(gè)強(qiáng)非線性問題，由于兩相界面的位置未知，因此界面上能量守恒條件也是非線性的．求解這類問題時(shí)，一般多用近似方法或數(shù)值方法求解[4-6]．?dāng)?shù)值方法求解時(shí)，通常應(yīng)用空間坐標(biāo)變換，把移動(dòng)區(qū)域模型轉(zhuǎn)化為固定區(qū)域模型進(jìn)行間接求解[7-8]．

本文以一維介質(zhì)的融化為模型，考慮在初始位置x=0處有一周期振蕩的熱源，融化從x=0處開始，固液交界面的移動(dòng)邊界X(t)沿x正方向移動(dòng)．令固態(tài)區(qū)域總保持恒溫并且恰好為相變溫度．而液態(tài)區(qū)域，通過空間坐標(biāo)變換，將移動(dòng)邊界轉(zhuǎn)化為固定邊界，采用顯式和隱式兩種有限差分格式數(shù)值模擬在融化過程中移動(dòng)邊界的運(yùn)動(dòng)和液態(tài)介質(zhì)內(nèi)溫度場(chǎng)的分布，進(jìn)而比較兩種有限差分方法的數(shù)值穩(wěn)定性、計(jì)算量和計(jì)算效率，為該類型的相變傳熱問題提供參考．

1 數(shù)學(xué)模型

令x=x′/X′(t′)，一維移動(dòng)邊界問題能夠通過空間坐標(biāo)的簡(jiǎn)單拉伸轉(zhuǎn)化為固定邊界問題[7]．這里x′是空間變量,0≤x′≤X′(t′);X′(t′)是移動(dòng)邊界位置;x是無量綱化的空間變量,0≤x≤1．

引進(jìn)無量綱變量和參數(shù)，在x=0處溫度隨時(shí)間周期振蕩的一維融化問題，其固定區(qū)域無量綱化模型可表示為[8]：

(1)

(2)

邊界和初始條件為

T(x=1,t)=0,t>0,

(3)

T(x=0,t)=1+εsin(ωt),t≥0,

(4)

X=0,t≤0,

(5)

2 兩種有限差分格式的建立

對(duì)時(shí)間進(jìn)行離散時(shí)，先取定一個(gè)小的時(shí)間t0作為計(jì)算的初始時(shí)刻，令tm=t0+mΔt，其中Δt為時(shí)間步長(zhǎng)，m=0,1,2,…．對(duì)空間區(qū)域x∈[0,1]進(jìn)行N等分，令xj=jh，j=0,1,2,…,N，其中h為空間步長(zhǎng)，并且x0=0，xN=1．

2.1 顯式有限差分格式的建立

微分方程組(1)-(5)的顯式有限差分格式為[8]：

j=1,2,…,N-1,

(6)

m=0,1,2,…,

(7)

(8)

(9)

X0=X(t0).

(10)

其截?cái)嗾`差為O(Δt+h2)．應(yīng)用凍結(jié)系數(shù)和Fourier穩(wěn)定性判定方法，上述顯格式在滿足如下條件時(shí)是數(shù)值穩(wěn)定的[9]：

2.2 隱式有限差分格式的建立

微分方程組(1)-(5)的隱式有限差分格式為：

j=1,2,…,N-1,

(11)

m=0,1,2,…,

(12)

(13)

(14)

X0=X(t0).

(15)

其截?cái)嗾`差為O(Δt+h2)．借助凍結(jié)系數(shù)的方法，可以證明該隱格式具有絕對(duì)的數(shù)值穩(wěn)定性[9]．

3 數(shù)值結(jié)果和討論

這里采用顯式和隱式兩種有限差分格式,對(duì)x=0處溫度周期振蕩(振蕩振幅ε=0.5，振蕩頻率ω=π/2)的條件下介質(zhì)融化過程中移動(dòng)邊界的運(yùn)動(dòng)和液態(tài)介質(zhì)內(nèi)溫度場(chǎng)的分布進(jìn)行數(shù)值模擬計(jì)算，并對(duì)兩種有限差分格式的數(shù)值穩(wěn)定性、計(jì)算量和計(jì)算效率進(jìn)行比較．這里所進(jìn)行的數(shù)值實(shí)驗(yàn)是用MATLAB 7.1編程、在系統(tǒng)為Microsoft Windows XP，CPU為Pentium(R)4/2.80 GHz，內(nèi)存為1 G的個(gè)人計(jì)算機(jī)上進(jìn)行的．

考慮純物質(zhì)的融化問題，對(duì)于熱源為定常溫度T(x=0,t)=1 (t>0)(即ε=0)的這種情況，微分方程組(1)-(5)的精確解為[10]：

(16)

(17)

其中λ滿足下面的超越方程

(18)

對(duì)于熱源周期振蕩的一維融化問題，需要說明的是為了避免t=0時(shí)的特殊情況，取t=0之后的某個(gè)時(shí)刻t0作為計(jì)算過程的初始時(shí)刻，應(yīng)用式(16)和(17)近似地對(duì)時(shí)刻t0處的溫度場(chǎng)和移動(dòng)邊界位置進(jìn)行初始化，這樣顯格式和隱格式兩種有限差分格式的數(shù)值計(jì)算過程就能順利向前推進(jìn)．

由顯式有限差分格式穩(wěn)定性的條件可知，顯格式對(duì)時(shí)間步長(zhǎng)的選取是有限制的且要求時(shí)間步長(zhǎng)非常小．這里取初始時(shí)刻為t0=0.01，空間步長(zhǎng)為h=1/N=0.1，振蕩振幅為ε=0.5.通過數(shù)值實(shí)驗(yàn)可知，當(dāng)時(shí)間步長(zhǎng)取Δt=0.000 02時(shí)，對(duì)于Stefan數(shù)分別取Ste=0.2、1.0和2.0，都可以保證顯式差分格式的數(shù)值穩(wěn)定性．若時(shí)間步長(zhǎng)取Δt=0.000 1時(shí)，對(duì)于Ste=2.0這種情況是符合上面所得顯格式穩(wěn)定性的充分條件，而對(duì)于Ste=1.0這種情況盡管不符合上面所給顯格式穩(wěn)定性的充分條件，但通過數(shù)值計(jì)算發(fā)現(xiàn)還是可以保證顯格式的數(shù)值穩(wěn)定性．這說明了上述所給顯格式的穩(wěn)定性條件是充分的，沒有包括時(shí)間步長(zhǎng)Δt的所有可能取值，主要是由于穩(wěn)定性分析時(shí)采用了“凍結(jié)系數(shù)”進(jìn)行預(yù)處理的結(jié)果．而隱式有限差分格式是一種絕對(duì)穩(wěn)定的差分格式[9]，對(duì)時(shí)間步長(zhǎng)和空間步長(zhǎng)的選取是沒有限制的．圖1顯示了在Ste=1.0、時(shí)間步長(zhǎng)Δt分別取0.1、0.01和0.001時(shí)應(yīng)用隱式有限差分格式所獲得的移動(dòng)邊界位置隨時(shí)間變化的關(guān)系曲線．可以看出所得的三條曲線近似程度非常高.特別地，當(dāng)時(shí)間步長(zhǎng)為0.01和0.001時(shí)，所對(duì)應(yīng)的兩條曲線完全重合；另一方面從數(shù)值上也進(jìn)一步表明了所構(gòu)造的隱式有限差分格式的絕對(duì)穩(wěn)定性．

圖2數(shù)值模擬了顯式差分格式和隱式差分格式對(duì)三個(gè)不同的Stefan數(shù)(Ste=0.2、1.0和2.0)所獲得的融化介質(zhì)移動(dòng)邊界的位置X(t)隨時(shí)間t的變化關(guān)系．從圖2(a)和圖2(b)可以看出，隨著Stefan數(shù)的增加，移動(dòng)邊界的運(yùn)動(dòng)速度也隨著增大，而且在x=0處周期振蕩的熱源對(duì)介質(zhì)相變速度的影響會(huì)在相變初期較短的時(shí)間段內(nèi)比較明顯．隨著時(shí)間的增加，在一段時(shí)間之后，移動(dòng)邊界的位置變化逐漸與時(shí)間的平方根成正比．此外，從圖2中還可以發(fā)現(xiàn)，使用這兩種差分格式所得的數(shù)值結(jié)果吻合地非常好，都能夠很好地反映在熱源周期振蕩條件下移動(dòng)邊界的位置隨時(shí)間的變化規(guī)律．

圖1 隱格式下不同時(shí)間步長(zhǎng)所對(duì)應(yīng)的移動(dòng)邊界的運(yùn)動(dòng)

圖2 不同Stefan數(shù)下兩種差分方法得到的移動(dòng)邊界的位置與時(shí)間的關(guān)系

圖3數(shù)值模擬了應(yīng)用兩種差分格式所得到的在不同時(shí)刻液態(tài)介質(zhì)內(nèi)的溫度分布．從圖3中可以看出，在x=0處周期振蕩的溫度強(qiáng)烈影響著液態(tài)介質(zhì)內(nèi)的溫度分布．通過比較圖3(a)和圖3(b)可以看出，這兩種差分格式所模擬的介質(zhì)內(nèi)溫度場(chǎng)的分布高度吻合．結(jié)合圖2和圖3，可推斷出采用這兩種有限差分格式對(duì)熱源周期振蕩條件下融化問題進(jìn)行數(shù)值計(jì)算都是可行的．

圖3 兩種差分方法得到的液態(tài)介質(zhì)內(nèi)不同時(shí)刻的溫度分布

下面從計(jì)算量和計(jì)算效率上對(duì)這兩種有限差分格式進(jìn)行比較．表1顯示的是h=0.1和Ste=1.0時(shí)兩種有限差分格式在計(jì)算到t=25.0時(shí)所需的計(jì)算量和計(jì)算時(shí)間．可以看出，對(duì)于顯格式，如果采用時(shí)間步長(zhǎng)為Δt=0.000 1，則計(jì)算到t=25.0這一時(shí)刻需要循環(huán)計(jì)算249 900次，計(jì)算量非常大；如果取時(shí)間步長(zhǎng)為Δt=0.001時(shí)進(jìn)行計(jì)算，數(shù)值結(jié)果不再收斂，出現(xiàn)了不穩(wěn)定的情況．然而對(duì)于絕對(duì)穩(wěn)定的隱格式來說，可以采用比較大的時(shí)間步長(zhǎng)．如果采用時(shí)間步長(zhǎng)為Δt=0.01，計(jì)算到t=25.0這一時(shí)刻只需要循環(huán)計(jì)算2 499次；若采用時(shí)間步長(zhǎng)為Δt=0.1進(jìn)行計(jì)算，仍然可以保證數(shù)值結(jié)果的穩(wěn)定性和收斂性，這時(shí)只需要循環(huán)計(jì)算250次，大大減少了計(jì)算量．此外，從表1中還可以看出，僅僅計(jì)算到t=25.0這一時(shí)刻，顯格式所需的運(yùn)行時(shí)間竟然達(dá)到13 303.137 805 s.而隱式差分格式在保證穩(wěn)定性和收斂性的前提下，時(shí)間步長(zhǎng)可以取Δt=0.1進(jìn)行計(jì)算，雖然每一次的循環(huán)計(jì)算都需要反復(fù)迭代才能獲得滿足誤差要求的結(jié)果，但只需5.242 976 s即可獲得t=25.0這一時(shí)刻的結(jié)果．因此，相比于顯式差分格式，隱式差分格式在介質(zhì)融化的理論分析中可以選取比較大的時(shí)間步長(zhǎng)進(jìn)行計(jì)算，能夠大大減少計(jì)算量和計(jì)算時(shí)間而更具有應(yīng)用價(jià)值．

表1 h=0.1和Ste=1.0時(shí)兩種差分格式計(jì)算量和計(jì)算時(shí)間的對(duì)比

4 結(jié)論

通過空間坐標(biāo)變換，把移動(dòng)區(qū)域模型轉(zhuǎn)化為固定區(qū)域模型，分別采用顯式有限差分格式和隱式有限差分格式求解了熱源周期振蕩條件下的一維融化問題．通過對(duì)融化過程中移動(dòng)邊界的運(yùn)動(dòng)及液態(tài)介質(zhì)內(nèi)溫度場(chǎng)分布的數(shù)值模擬發(fā)現(xiàn)這兩種差分格式的數(shù)值結(jié)果高度吻合，但是在數(shù)值穩(wěn)定性、計(jì)算量和計(jì)算效率方面隱式差分格式明顯優(yōu)于顯式差分格式．