趙妍，令鋒

(1.內(nèi)蒙古工業(yè)大學(xué)理學(xué)院，內(nèi)蒙古呼和浩特010051；2.肇慶學(xué)院計算機(jī)學(xué)院，廣東肇慶526061)

熔化凝固問題的Keller盒式格式

趙妍1,2，令鋒2

(1.內(nèi)蒙古工業(yè)大學(xué)理學(xué)院，內(nèi)蒙古呼和浩特010051；2.肇慶學(xué)院計算機(jī)學(xué)院，廣東肇慶526061)

建立了熱參數(shù)依賴于溫度的熔化凝固問題的顯熱容模型的Keller盒式格式，運用凍結(jié)系數(shù)法證明了該格式的無條件穩(wěn)定性，并給出了數(shù)值算例.

顯熱容法；盒式格式；穩(wěn)定性

數(shù)值求解熔化凝固問題的顯熱容模型[1]的控制方程具有強(qiáng)非線性性，其熱參數(shù)是依賴于溫度的函數(shù).在數(shù)值求解中，發(fā)生相變的溫度區(qū)間的選取對計算結(jié)果有較大影響.為克服這一缺點，1985年，Hsiao對顯熱容模型做了改進(jìn)[2].Keller盒式格式是Keller在1970年提出的一種差分格式[3].1984年，Meek和Norbury用Keller盒式格式求解了非線性移動邊界問題[4]；2008年，Hamzah等給出了Keller盒式格式的并行算法[5].本文中，筆者在前人工作的基礎(chǔ)上，對熔化凝固問題的顯熱容模型做進(jìn)一步改進(jìn)，用Keller盒式格式離散其控制方程，并分析該格式的無條件穩(wěn)定性,最后給出數(shù)值算例.

1 顯熱容模型

根據(jù)顯熱容法，熔化凝固問題的控制方程可表示為

初始條件為

邊界條件為

式中：

C(T)是依賴于溫度的有界間斷函數(shù)，Tm是相變臨界溫度.

2 Keller盒式格式

故Keller盒式格式的截斷誤差為R＝O(τ2＋h2).

文獻(xiàn)[2]657－661的研究結(jié)果表明，直接應(yīng)用式(4)會導(dǎo)致數(shù)值誤差或解的不穩(wěn)定性，為此，筆者依據(jù)文獻(xiàn)[2]657－661的方法對熱容量做進(jìn)一步改進(jìn).由文獻(xiàn)[2]657－661(見圖1)知，在(i，j)點的熱容量C(Ti，j)用(i－1，j)，(i＋1，j)，(i，j－1)，(i，j＋1)4個點熱容和的平均值來代替,即

依據(jù)同樣的原理(見圖2)，在(j－1/2,n－1/2)點的熱容量C(Tnj－－11//22)用(j－1,n)，(j,n)，(j－1,n－1)，(j,n－1)4個點熱容和的平均值來代替,即

其中：Tn－1/2表示節(jié)點(j－1/2,n－1/2)處的溫度，C(Tn－1/2，Tn)表示修正熱容量.

圖1 改進(jìn)的顯熱容法的節(jié)點示意圖

圖2 Keller盒式格式的節(jié)點示意圖

3 Keller盒式格式的穩(wěn)定性分析

定理1差分格式(12)是無條件穩(wěn)定的.

證用“凍結(jié)系數(shù)”法[6]討論其穩(wěn)定性.

使用記號(7)，將式(12)展開得

其中：λ＝τ/h2.

令Tnj＝vneikjh，代入式(15)得

將式(16)整理后得

4 數(shù)值算例

為驗證本文算法的有效性，以下用Keller盒式格式求解存在解析解問題[7](18)～(21)的數(shù)值解：

其解析解為

由積分公式及式(18)與(19)，式(21)中的第2式可以寫為如下形式：

控制方程(25)定義在固定邊界0＜σ＜1,0＜τ＜H上.

邊界條件(19)變形為

初始條件(20)變形為

邊界條件(21)中的第1式變形為

式(23)變形為

利用Keller盒式格式，控制方程(25)可以寫為如下形式：

控制方程(30)離散為下列2個差分方程：

在σ＝0處的邊界條件可近似表達(dá)為

初始條件(27)離散為

邊界條件(28)離散為

對兩相界面的守恒條件（式(29)）應(yīng)用矩形公式，得

求解方程組(31)～(36)時，取Δσ＝0.01，Δτ＝0.05，N＝100，H＝5，J＝100.

迭代過程的一般步驟如下：

1)先給出試探的s(τn＋1)；

2)用TDMA算法求解方程(31)～(35)，求得vnj＋1，unj＋1(j＝0,1,2,…,J－1)；

3)根據(jù)求得的unj＋1，由式(36)計算出新的s(τn＋1)；

4)用新的s(τn＋1)重復(fù)步驟2)，直到滿足收斂條件

圖3u(3,t)時數(shù)值解和精確解之差

圖3 為點(3,t)的數(shù)值解和精確解之差.從圖3中可以看出，本文得到的數(shù)值解非常接近于精確解，表明筆者提出的方法是求解熔化凝固問題的有效方法.

[1]BONACINA C,COMINI G,FASANO A,et al.Numerical solution of phase-change problems[J].Heat Mass Transfer,1973,16：1 825-1 832.

[2]HSIAO J S.An efficient algorithm for finite-difference analyses of heat transfer with melting and solidification[J].Numerical Heat Transfer,1985(8)：653-666.

[3]KELLER H B.A new difference scheme for parabolic problems,in Numerical Solutions of Partial Differential Equations[M]. 2nd ed.New York：Academic Press,1971：327-350.

[4]MEEK P C,NORBURY J.Nonlinear moving boundary problems and a Keller box scheme[J].SIAM J Numer Anal,1984,21：883-893.

[5]HAMZAH N,ALIAS N,AMIN N S.The parallelization of the Keller box method on heterogeneous cluster of workstations[J]. Journal of Fundamental Sciences,2008(4)：253-259.

[6]陸金甫,關(guān)治.偏微分方程數(shù)值解法[M].2版.北京：清華大學(xué)出版社,2004：57-59.

[7]FURZELAND R M.A comparative study of numerical methods for moving boundary problems[J].J Inst Maths Applics,1980, 26：411-429.

Keller Box Scheme for Melting and Solidification Problem

ZHAO Yan1,2,LING Feng2
(1.College of Sciences,Inter Mongolia University of Technology,Hohhot Inter Mongolia010051,China; 2.School of Computer Science,Zhaoqing University,Zhaoqing Guangdong526061,China)

The Keller box scheme for the mathematical model of melting and solidification problem with temperature-dependent heat parameters is established.The unconditional stability of the scheme is analyzed by using the freeze coefficient method.And a numerical example is presented.

apparent heat capacity method;Keller box scheme;stability

O241.82

A

1009－8445（2010）02－0001－05

(責(zé)任編輯：陳靜)

2009－10－09

廣東省自然科學(xué)基金資助項目(04011600)

趙妍(1983－),女,黑龍江哈爾濱人,內(nèi)蒙古工業(yè)大學(xué)與肇慶學(xué)院聯(lián)合培養(yǎng)碩士研究生.