曲良輝 邢 琳 錢德亮 令 鋒

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

曲良輝 邢 琳 錢德亮 令 鋒

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

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

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

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

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

滿足的條件為

引進(jìn)無量綱量

定空間步長方法

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

因此，應(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)用定空間步長方法，圖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)溫度的變化曲線，其中空間變量步長固定為Δr=0.001。表1對定空間步長方法所得到的不同時(shí)刻處移動(dòng)邊界的位置與文獻(xiàn)中的相應(yīng)結(jié)果進(jìn)行了比較。從表中可以看出，在Stefan數(shù)Ste=0.2和邊界條件f（ t）=0的情況下，定空間步長方法所得的結(jié)果與文獻(xiàn)中的焓法和攝動(dòng)法所得到的相應(yīng)結(jié)果非常接近，其絕對誤差均不超過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)用定空間步長方法所得到的移動(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)溫度的變化曲線，其中使用的空間變量步長為Δr=0.001。從圖3（a）中可以看出，在一段時(shí)間之后，移動(dòng)邊界的位置隨時(shí)間的變化幾乎保持線性的運(yùn)動(dòng)趨勢，也就是說移動(dòng)邊界的運(yùn)動(dòng)速度幾乎不變。從圖3（b）中可以看出，相變區(qū)域內(nèi)的溫度場在快速變化，當(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é)束語

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

DOI：10.3969/j.issn.1001-8972.2015.15.005