曲良輝 邢 琳 錢德亮 令 鋒

低溫邊界條件下柱型凝固問題的數(shù)值研究

曲良輝 邢 琳 錢德亮 令 鋒

通過考慮移動(dòng)邊界每次向外推進(jìn)固定距離時(shí)所需要的時(shí)間，構(gòu)造定空間步長(zhǎng)方法求解了兩種不同低溫邊界條件下的向外柱型凝固問題，數(shù)值模擬了相變過程中移動(dòng)邊界的運(yùn)動(dòng)及介質(zhì)內(nèi)溫度場(chǎng)的變化。數(shù)值結(jié)果表明采用定空間步長(zhǎng)方法求解柱型凝固問題是數(shù)值穩(wěn)定的，并且具有較高的精度，對(duì)于Stefan問題的研究具有重要的參考價(jià)值。

Stefan問題也稱為移動(dòng)邊界問題，其特點(diǎn)是區(qū)域內(nèi)存在著一個(gè)隨時(shí)間運(yùn)動(dòng)的邊界，如冰的融化、金屬鑄件的凝固、微粒擴(kuò)散和種群繁衍等。求解Stefan問題，由于問題的非線性性質(zhì)和移動(dòng)邊界的未知位置，只有一些特殊情況才能獲得解析解，通常只能用近似方法或數(shù)值方法求解。對(duì)于柱狀模型下的相變傳熱問題，1981年Lunardini借助熱平衡積分技巧，應(yīng)用等效熱擴(kuò)散率方法獲得了具有定常熱源邊界條件的柱型相變問題的一個(gè)近似解析解。2009年，Caldwell等應(yīng)用焓法、邊界固定法、攝動(dòng)法和熱平衡積分法分別求解了柱型相變傳熱問題，并就定常邊界條件的情況，對(duì)數(shù)值結(jié)果進(jìn)行了分析和比較。

本文考慮在低溫邊界條件下，通過計(jì)算移動(dòng)邊界每次向外推進(jìn)固定距離時(shí)所需要的時(shí)間，構(gòu)造了一種定空間步長(zhǎng)方法數(shù)值求解向外柱型凝固問題，并對(duì)移動(dòng)邊界的運(yùn)動(dòng)和固態(tài)介質(zhì)內(nèi)的溫度場(chǎng)進(jìn)行了數(shù)值模擬，分析了低溫邊界條件對(duì)相變區(qū)域內(nèi)的溫度變化及移動(dòng)邊界運(yùn)動(dòng)的影響。

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

考慮初始溫度恰好為相變溫度Tf的某液態(tài)介質(zhì)由于邊界處溫度的變化所導(dǎo)致的向外柱型凝固問題（圖1），在相變過程中體積與物性均不發(fā)生變化的情況下，該問題的數(shù)學(xué)模型為

滿足的條件為

引進(jìn)無(wú)量綱量

定空間步長(zhǎng)方法

記Δr 為固定的空間網(wǎng)格步長(zhǎng)，tj（t0=0）是對(duì)應(yīng)于移動(dòng)邊界位于Rj=1+jΔ r （j=0，1，…）處的時(shí)間?？紤]相變界面移動(dòng)到RN=1+NΔ r 處所對(duì)應(yīng)的時(shí)間tN，此時(shí)各節(jié)點(diǎn)rj=1+jΔ r 處的固相介質(zhì)的溫度記為由于控制方程可以簡(jiǎn)化為

因此，應(yīng)用向后差分格式離散時(shí)間導(dǎo)數(shù)項(xiàng)和中心差分格式離散空間導(dǎo)數(shù)項(xiàng)，則控制方程的離散形式為

Stefan條件（7）的離散形式為

引入兩個(gè)變量a=1Δr 和b=ΔrΔtN-1，式（11）和（12）分別變?yōu)?/p>

邊界條件（8）-（10）的離散形式分別為

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

圖1　向外柱型凝固示意圖

在Stefan數(shù)Ste=0.2和f（ t）=0的情況下，應(yīng)用定空間步長(zhǎng)方法，圖2（a）數(shù)值模擬了移動(dòng)邊界位置隨時(shí)間變化的運(yùn)動(dòng)曲線和圖2（b）數(shù)值模擬了當(dāng)移動(dòng)邊界分別向外移動(dòng)到R（ t）=1.2，1.5和1.8時(shí)，相變區(qū)域內(nèi)溫度的變化曲線，其中空間變量步長(zhǎng)固定為Δr=0.001。表1對(duì)定空間步長(zhǎng)方法所得到的不同時(shí)刻處移動(dòng)邊界的位置與文獻(xiàn)中的相應(yīng)結(jié)果進(jìn)行了比較。從表中可以看出，在Stefan數(shù)Ste=0.2和邊界條件f（ t）=0的情況下，定空間步長(zhǎng)方法所得的結(jié)果與文獻(xiàn)中的焓法和攝動(dòng)法所得到的相應(yīng)結(jié)果非常接近，其絕對(duì)誤差均不超過0.001，從而表明了該數(shù)值方法求解柱型凝固問題的可行性和有效性。

圖2　在Ste=0.2和f（ t）=0的情況下移動(dòng)邊界的運(yùn)動(dòng)和相變區(qū)域內(nèi)的溫度分布

表1　在Ste=0.2和f（ t）=0的情況下移動(dòng)邊界位置的比較

在Stefan數(shù)Ste=0.2和f（ t）=-t2的情況下，圖3（a）數(shù)值模擬了應(yīng)用定空間步長(zhǎng)方法所得到的移動(dòng)邊界位置隨時(shí)間變化的運(yùn)動(dòng)曲線，圖3（b）數(shù)值模擬了當(dāng)移動(dòng)邊界分別向外移動(dòng)到R（ t ）1.2，1.6和2.0時(shí)，所得到的相變區(qū)域內(nèi)溫度的變化曲線，其中使用的空間變量步長(zhǎng)為Δr=0.001。從圖3（a）中可以看出，在一段時(shí)間之后，移動(dòng)邊界的位置隨時(shí)間的變化幾乎保持線性的運(yùn)動(dòng)趨勢(shì)，也就是說移動(dòng)邊界的運(yùn)動(dòng)速度幾乎不變。從圖3（b）中可以看出，相變區(qū)域內(nèi)的溫度場(chǎng)在快速變化，當(dāng)移動(dòng)邊界的位置從R（ t）=1.6向外推進(jìn)到R（ t）=2.0時(shí)，相變區(qū)域內(nèi)的溫差從大約1.8增加到了約4.3，這是受低溫邊界條件影響的結(jié)果。

圖3　在Ste=0.2和f（ t）=-t2的情況下移動(dòng)邊界的運(yùn)動(dòng)和相變區(qū)域內(nèi)的溫度分布

結(jié)束語(yǔ)

通過考慮相變界面每次向外推進(jìn)固定距離時(shí)所需要的時(shí)間，構(gòu)造了定空間步長(zhǎng)方法用來(lái)求解向外柱型凝固問題，并模擬了相變過程中移動(dòng)邊界的運(yùn)動(dòng)和固態(tài)介質(zhì)內(nèi)溫度場(chǎng)的變化。通過與已有結(jié)果比較可知，應(yīng)用定空間步長(zhǎng)方法求解向外凝固問題是可行的，并且具有較高的精度，從而表明了定空間步長(zhǎng)方法在求解Stefan問題時(shí)所具有的重要應(yīng)用價(jià)值。

DOI：10.3969/j.issn.1001-8972.2015.15.005