Stefan问题的数值解法

‘岌犷花歹阅乒飞歹花步~认川噶司瞥进展 St ef an 问 题 的 数 值 解 法 忻 孝 康 (复旦大学应用力学系 ) 引 言 在自然界和各种工业生产中 , 大量存在着具有相变的热传导问题 , 例如 : 地球两极的极地冰 块的熔解(M elt in g )和固化 (S o lid ifie a t io n ) , 炼钢中铁块熔化 及浇铸时铸件的固结 (凝结 ) ; 熔岩 的冷却 、 固结问题 , 飞机穿过云层时机翼前缘的结冰问题以及制冷工业中冰的熔解和固结问题 , 都是广泛存在着的这类问题 . 这类具有相变 (从液相变到固相或从固相变到液相 )的传热问题是传 热问题中一个特殊的方面 。 这类问题的特点是 : 控制方程是热传导方程 , 而区域内含有一个或 。 个运动 (移动)的分界面 (Int e r fac e ) , 在该分界面上放出或吸收热量 . 这 类问题统称 为 S te fan . 问 题 , 以纪念上世纪末曾对此作过系统研究的科学家 J . Ste fa n( G . L a m e 和 B . P . Cl a pe yr o n 最 早于 188 1一 1 8 8 3 年曾研究过半无 限大液体区域(x > 0) , 由于 二 “0 边界上的低温 引起液 体 凝 固 问题 , 他们得到过固结厚度 S (t )正比于训了的结果 , 但没有找到比例 常 数 . 188 9一18 9 1 年间 , J . S te fa n 研究了极地冰块的熔解和土壤冻结问题 , 给出了级数近似解和相似解 S (t) = ZA 了丁 , 特别是给出了分界面上的热平衡条件 : 、.了、J1 声夕矛.瓜,‘O “. , _ 0 “: _ 、 _ d 名(t) 几 . 一了二丫 一 凡名冲I丫了 一 几f, 一j 玉一 .口劣 口劣 “ 石 其中二 = : (t )是分界面位置 , 入为单位质量的熔解潜热 , ‘ , : 分别为固相 、 液相的温度分布 , 从 , 别为固 、 液相的热导系数 t又‘“1) 。 以后 , 人们发现 , 在许多扩散问题中也存在类似的现象 , 其控制方程为扩散方程 C , = (D C , )。 , 也存在一个运动的未知边界 , 在该边界上成立两个边界条件 (一个边界条件类似于(1) 式 ) . 从而 , 把此类扩散问题也统称为 st efa n 问题 . 因此 , 目前人们都采用如下的定义 [“] : 经典 St e fa n 问题 : 控制方程是扩散 (或扩散对流)方程的某个初边值问题 , 区域内含有 一 个 或几个运动边界 (在该边界上满足两个边界条件)的问题 , 统称为经典 S te fa n 问题 . 广义 St e fa n 问题 : 控制方程是一般的抛物型方程 (例如 , “。= au 二二 + 加二 + 阅 + d) 的 某个初 - 边值间题 , 区域内含有一个或多个运动边界的问题 , 就统称为广义 St efa n 问题或 S t ef a 攀 型间题 (S te fa n T vPe P r o b le m ) 其它如多孔介质中液体的蒸发 、 化学溶液里反应物的沉淀〔‘1肌肉中氧气的扩散 [ “]等类问题 , 均可归结为 S te fa n 问题 . St ef a n 问题 数学 数学高考答题卡模板高考数学答题卡模板三年级数学混合运算测试卷数学作业设计案例新人教版八年级上数学教学计划 上是一个强非线性问题 , 即使其控制方程是线性 的 (例如线性热传导方程 ) , 但是由于内部或区域一部分边界是未知的 (运动的 ) , 其边界条件是非线性的 (见下节 ) , 因此 , 要 寻求解析解是相当困难的 , 特别是二维以上的高维问题 . 文献〔习给出了一维问题各种情况下的详 一 89 一 细的解析解结果 . 虽然, 直至 1 9 79 年还有人对一维头e fa 。 问题的辐射边界条件情况给出了级数 形式的解析解 (J . Fl ui d M ec h a ni cs ), 但是总的趋势是倾向于数值解法和近似解法 . 对 St ef a n 问题的理论研究 , 在数学上也十分有趣 , (广义 )S te fa n 问题的存在性和唯一性问题 , 至今仍吸引着许多数学家在进行工作 , 这方面的工作可参见文献〔6一8〕. 除了严格的解析结果之 外 , 许多近似 (半解析 )方法也是相当重要的 , 可借的是许多近似方法在高维问题的推广比较困难 ; 即使可行 , 要得到结果仍需数值计算。 这些近似和半解析的方法有 : 积分方法 [。l(这有点象 边 界 层理论中的 K a r m a 卜P o hlh a u s e n 方法 )、 变分方法 [’“} 、 G r e e n 函数方法I“ ] (这与边界积分方 程 方法一样 )以及摄动方法t‘“]等 . 由于篇幅关系 , 这里我们不准备对这些近似或半解析方法进行 介 绍和分析 , 而把讨论的重点放在纯数值方法上 . (应该指出的是 : G ree n 函数方法或运动的 源 汇 分布方法 , 实质上是边界积分方程方法 , 因此同样可以用所谓 ‘边界元 ” 方法 [ ‘“〕来数值求解 , 故 也可算是一种数值方法 ) ’ 5 te fa n 问题的数值解法 , 主要有: 有限差分法t“ , ‘”] 、等温边界移动方法 (T h e r m al M ig rat io n M e th o d ) [ ‘e , ‘7 1 、 贴体坐标系方法 (F ittfn g B o d y C o o r d in a te Sy ste m M e th o d ) [且. , ’. ] 热容量或 热烩方法 [ 2 “ , “气 运动网格的有限元法I“’〕变分不等式方法 (法国 )及变区域上的变分有限元法t Z“, “月 等 . 由于分界面位置是未知的 , 边界条件是非线性的 , 因此到 目前为止尚没有一个所谓 ‘最好 ” 的方法 。 往往某些方法对某些问题特别有效 ; 有些方法对维数的限制较少 , 有的则限制较多 . 总 之 , 各有各的优缺点 , 要看问题的特性如柯 。 本文的 目的是介绍 目前国内外数值求解 St e fa n 问题的各种主要的数值方法 , 以供初步 涉 及 St ef a n 问题数值研究的同志有一个大致的了解 . 二 、 St ef a n 间题的数学描述 为了便于下面的介绍和叙述 , 让我们以一维杆 (棒)的熔解或凝固(固结或冻结)问题为例来看 St e fa n 问题的数学提法 . TTI爪几 考虑一根半无限长 (0 《劣< co ) 的金属棒 , 初 始温度为均匀温度 T 。, 设金属的熔解温度为 T脚 假设 T 。< T⋯ 如果突然在 棒端 二 = O 处加热 至温度 丁: (> T砂 且一直保持此温度 ; 于是 , 该棒从左 端(二 = 0) 开始逐渐熔解 , 形成了一个从左向右移动 的熔解分界面 : 二 二 S (t ) 。 它分割着左面的液相区 域和右边的固相区域 . 如果不考虑固相和液相的密 度差 , 即 p = p . = p‘(次后均用下标 s 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示固相 , 矛 表示液相 )(见图 1 ) . · - 那么 , 该问题的数学提法是 : 液相 武t) 固相 p。: ·争 ·斋卜会)一 之 。c . 鲁 一斋(xa会) 0 丈戈 < S (t ) (Z a ) S (t) < 戈 < 切 (Z b ) (其中 C 。 , ‘, 棍 , . 分别为固相和液相的比热和导热系数 ) 一 90 一 初始条件 : T . (二 , 0 ) = T 。 0 < 戈( oo (3 a ) S (0 ) = 0 (3 b ) 固定端边界条件 : T : (0 , t) = T : , t) 0 (4 a ) T 一 (+ co , t) = T 。, t) 0 (4 b ) 留下来的是运动分界面 劣 = S (f )上的边界条件 。 一个边界条件是很显然的 , 即温度连续条件 : T 一 (S (t ) , t) = T z(S (t) , t ) = T 喊 (5 ) 分界面上的另一个边界条件可以从分界面上的热平衡条件来得到 . 大家都知道 , 任何 物 质 从 固 相熔解成液相必须吸收热量 , 而由液相变成固相时又必须释放出热量 (与吸收的 热 量 相 等 ). 若 用 入表示单位质量物质的吸收 (或释放 )的熔解潜热 , 则单位时间内由固相熔解成液相所需要吸收 的热量应为 : 久p d s (t) d 矛 ‘ 、、 ~ 二 山 , 一 。 、 ‘二 , 止 , J T : 山 , 、 , 一一 ~ ~ 。 . , . 一一 。 、 , a T , ~ _ ” _早 1止浏 I川 曰 , 出限怕汗 / \汀习只卜更;刀 一 凡 l we 万一二一 ; 出万开圆 l叫回怕书出划恐夏刀 一 K. ~ 5 二 , 囚此 , 夕了口 劣 口工 界面上的热平衡条件为 : 一 * :孕 一(一 、, ,孕)· 、p擎L, 荡 、 口不 / “ 不 、.孕 一 、:孕 一 * p塑U 荡 U 笑 “ 不 在 二 = S (t)上 (6 ) 这就是 J . S te fa n 最早导出的分界面热量平衡关系式 (图 2 ) . 在高维的情况下 , 类似的有 : aT , , aT , k . 羊鑫弃一 无:岑牛 = 入p 。。’ 一‘ 0 ” ‘ ” a n ‘ ’尸一 (7 ) }寿力 一 91 一 其中 n 是由固相指向液相的分界面单位法向 , 。 , 是分界面移动的法向分速度。 这样 , 方程式 (2) 初始边界条件(3 一6) 就构成了上述熔解问题的数学提法 (对于固结问题几乎 有同样的提法 , 参见 〔1」、 〔2〕) . 这个问题通常称为一维两相(T w o Phas e) St e fa n 问题 。 作为这 个问题的特例 , 如果初始时刻棒处于均 匀的熔解温度 T 。下 , 即 T 。二 T 。 , 那么方程(Zb) 自动 消失 (即 T : 二T . ) , 这种问题称为一维单相 (O。一Ph a “ )S tef a n 问题 。 (注意 : (S te fa n 问题的 ‘相 ” 说法 , 不是由物理上涉及一相或两相来划分的 , 而是由控制方程的个数来决定的 . ) 类似地 , 我们可以得到高维熔解问题的数学提法 t“名] . 例如 , 一个固体物质初始时刻 为 均 匀 温度 T 。(< T动 , 初始占有区域为 D , 边界记为 aD 。 t = O 时 , 突然在外边界 aD 上加温 至 T ‘并 一直保持此温度 (几 > T 动 , 这样在D 内就产生了一个移动的分界面 r , 它分割着液相和固相区域 (图 3 ) . 该间题的数学提法是 : pc誉 一 v p c. 鲁 一 v (无. V T ‘) 在D ‘(t )内 (s a ) (k * V T , ) 在D . (t)内 (s b ) 初始条件 : T . (劣 , 0 ) = T 。 及 (o ) 边界条件 : T : (劣 , t) = T : 分界面 上条件 : 在D : (0 )内 = 0 (g a ) (g b ) 在aD 上 (10 ) D == D. + 肠 图 3 T : 】r 〔, ) = T . ! : “ 〕= T , . 、, 琴一 、‘擎 = 、p 。。口刀 口打 (1 1 ) (1 2 ) 条件(11) 的意义同 (7) 式。 这是一个高维 (二或三维)的两相 S te fa n 问题 . 类似地 , 当 T 。二 T . 时 , 问题化为单相 Ste fa n 问题 , 这时方程 (8 b) 不必考虑 . 在这一节的最后 , 让我们再来看一下 St e fa n 问题的非 线 性 。 例如 , 对 于一维 两 相 st e fa n 问题 (2) 一 (6) 。 乍看上去 , 方程式 (助 。 初始和边界条件 (3) 一(6) 似乎都是线性的 , 问题似 乎 应 该是线性的 。 实际上 , 由于分界面位置 x = S (f) 是未知 、 待求的 , 关系式 (6 )实质上是一个非 线 性关系式 . 由 (石)式知 : T ‘ (S (t) , t ) = T . T , (S (t ) , 卖) = T , } (1 3 ) 对此两式取全导数 , 则有 [鲁 d X + 誉 d ‘}二 一 : (, ) · }会 d 二 + 鲁d ‘}二 一 : (, ) 二 ” 一 9 2 一 戈 二 S (t) 一誉机”而 一一一兰生 } d t }戈 把它代入 ( 6) 式 , 则有 a T : { = 一 一互{ = S ( t ) 些卫{ d 戈 }二 == S ( t ) aT : aT , , aT , 、 a t 凡一玉, , 一 ‘ i , 二: 一 = 一 p 人 一可不苏一一 = U 工 口万 口 1 止 (14 ) a x 这时 , 非线性关系就变得更加清楚了 . 因此 , 分界面上关系式 (助实际上是一个非线性关系式 . 因 此 , S te fa n 问题实质上是一个非线性问题 , 这也就是为什么它至今仍吸引许多人致力于这个问题 研究的一个原因 。 三 、 Stef a n 问题的数值解法综述 由于 St e fa n 问题的非线性特性 , 解不能进行叠加 , 因此要寻求问题的精确解相当困难 . 对一 维的问题 (单相或两相 ) 已给出了许多封闭形式的解析解 I ’闪 , 对二维 问 题 , 大 约 只 有 B ol ey & Y ag od a 的解析工作t Z盯 . 因此 , 高维的 S t e fa n 问题多数都是用各种近似方法和 数 值 方 法 加 以 解决的 。 限于篇幅 , 这里只对现有的一些主要的数值方法作一介绍。 理论上说 , 一些离散化的方法 (有限差分法 、 有限元法 、 边界元法 、 加权残数法等)都可以直 接应用于 S t e fa n 问题的求解。 但由于运动边界的存在 , 有限元法和边界元法的应用 目前仍 相 当 少 [ 2‘] . 数值方法中仍以有限差分法为主 。 有的是直接在固定网格或变时间网格上应用差 分 格 式 〔’‘” S J ; 有的是作了自变量或因变量变换后应用各种差分格式 〔‘“一 ““J . 在有限元方法中 , 目前看 来 变区域上的变分有限元方法是一种较有前途和较精确的方法 . 现按照方法的类型 , 对一些方法加 以简述之 。 (一 )有限差分法 1 . 固定步长的方法 这 时 , △二 和 △t 是固定的 , 因此分界面位置 x 二 S (t )一般说来不会落在网格节点上 (图 4 ) 。 X = 名( t) 耘+ - ‘(‘ + t ) 汁止击扒( : . 》 文献 [ 1 4〕考虑T 如下 的一维 单相 S t e fa n 问题 : . , = ‘二 , 0 < 义 ( S ( t ) “二 ( 0 , t ) = 一 1 , t> 0 u ( S ( t ) , t ) 二 0 , t ) 0 ( 1 6 ) S ( 0 ) ‘ n u 二 ( S ( t ) , t ) = d s ( t ) d t ‘ 他们首先把 (1 5 )最后一式化成如下等价的 积分式 : s ( , 》一卜丁:“‘ · (二 , ‘, d 、 ( 1 6 ) 一 93 一 然后 , 取 对内点 、 r _ 「夕(t, ) 1 工V 石 . ,- 石一~ 一 l L 已戈 」 (取整数) = 1 , 牙二N 一 1 , 列出方程的隐式格式 : “梦小 ’一 。犷_ “犷李圣+ u 梦士1一 加梦+ ‘一五t 一一 △工‘ 一 = 1 , 2. ·N 一 1 (灯)由于 “(S (t。十 : ) , t。 , , ) = 0 , 隐式差分格式 , “井+ ‘一 “井 △t 左端边界条件可近似为: 对 (N 一 1 , N , S (t . , : ” 三点 , 利用抛物线擂值 , 可以得到N点上的 = -红了 _ 一丛兰 _ _ _ △戈 L S (t。十 x ) 一 (N 一 1 )△劣 “黔‘ S (t 。 + ; ) 一 N △x } (18 ) 衅 + ’一 心 备‘= 一 △岑 (1 9 ) 最后 , 对(1 助式利用矩形积分公式 : S (t 。 + : ) = f 。 十 , 一 兄 “梦‘ ’ △二 - 1 , 。 , . 百L。‘才。 , : ) 一 (N 一 1)△劣〕“井+ ‘ (20 ) 从而可以得到如下的迭代过程 : ¹ 先给定一个猜测值 5 . ( t . , :》; 由差分方程 (竹)一 (19 )解一个三对角方程组 , 解出 “扩‘⋯ “扩 ‘: 再按 (加 )式计算出新的 S ! ( t 。+ ! ) , ’ !- 重复步骤 º , 直至 !S , + , ( t。 十 ; ) 一 S伙r . + : ) I < 。 。 l 这个方法在一维问题时是可行的 , 但内点的精度 O (△二, ) 与万点 ‘” 精度 ( O ( △幻 )不一致 . 山 . 2 . 变时间步长 Af 的办法 . 由于内点与边界点 (附近 )精度不 ‘ 一致 , 因此 D o u g la s一all ie[ “,设计了一种△x 固定而改变 △厂的 办法 . 选取步长 △t . 使得 S (t . + : ) 一 S ( t. ) 正好为一个 △x 的距离 ; 这样每走一个时间 △人 , 点数正好增加一个 (图 6 ) , 而离散方程 组为 : N 一 I N º¼» 扮 + 1 0 I N 一 1 / N 些士址 .竺兰= 川打 + “_出 一加 , + ‘△ t △另2 = 1 , 2. 二 N ”空令 ‘一 “忿十 ‘= 一 △x u哭李圣二 0 (21 ) S ( t . , : ) 二 t。 + △f。 一 艺 “梦“△二 = s ( t . ) + △x (22 ) 从而迭代可如下进行 : ¹ 选迭定一个 △九 , 由(21) 解出 : 李+ ‘; º 由 (22 )式算出一个新的△t 。 ; » 迭代到 I△代+ ’一 △坟 !< “ . 3 . 高维问题的差分法 [ ‘“I 上述固定空 间网格的办法 , L az ar id is 已发展到高维问题其基本思想类似干〔14〕, 例如对 内 点 (远离分界面 ) , 可以用通常的五点格式离散方程 (8 a) ; 一 9 4 一 皿塑二乙达一 二 。 , . (凸t , ’J 、T 梦东圣, , + T r士圣, , 一 ZT 梦亨笠△x 艺 沪呼+ 飞 . _ Q , 考1 小 7’ 呼+ 芝 . 、+ 立匕二竺一公少一曰兰) (2 3 ) 其中 。 =今为热扩散系数 , 设为常数 .尸七 对于接近分界面上的点 (如图 6 中的尸点 ) , 可以用 如下的时间后差离散公式 : 勺, . + 二_ 口, .占 , 盛 p _ ~一一 — 一 飞J . . .△ t r (豁 +等), (2‘, 匆 为了计算 (V Z T ), 值 , 可以采用尸点的双二次 (抛物面 ) 插值公式 : p l气 之肠. T (戈 , 夕) = T , + A (, 一 戈一) + B (v 一 , , ) + C (戈 一 戈, ) + D (夕一 梦, ) + E (戈 一 劣 , ) (g 一 v , ) 利用 T 尸 , T 。, T * , T , 及 T 二 和 △戈 , △y , 乙x (2石) 各夕可 图 6 以完全确定 A 一E 的值 。 从而 (V ZT ), = ZC + ZD (2 6 ) 这样差分格式完全确定 . 为了达到较好的精度 , 作者采用了很复杂的判别过程来计算整个问题 。 可以看到 , 这个方法的缺点在于分界面不位于 网格节点上 , 因此需要一个很复杂的插值过程 , 程序相当繁复 。 (二 )变换成固定 区域 的方法 利用 自变量或因变量变换 , 或 自变量与因变量位置的互换 , 对某些 St o fa n 问题的求解 是有 利的 , 甚至可使变区域的间题化为 固定区域的问题 。 目前大致有这么几类方法 : 1 . 等温边界移动方法 (T h e r m a l M ig r a tio n M e th o d ). 这主要是 C ra n k 等人的工作I’“, 川 。 譬如考虑如下的一维单相 S tef a n 问题 : “, = “二二 , 0 < 戈 < S (f) 。 } = “。 (> 0 ) S (0 ) = 0 。 ! 。。, , = 0 a u Ox {一。‘, , (2 7 ) = _ 日止丝 d t 注意到温度在边界上达极值的现象 (数学上可以严格证明) 。 分界面是一条等温 ( = 0) 线 。通常我们 都理解为 “ = : (x , t) , 如果我们把 “ 看作是自变量 , x 看作是因变!量 , 即 二 = 以“ , t) , 那 么 有 什么新的发现呢? 这时 , 很有意思的是 : 在 (。 , t) 平面上区域是一个固定的半无限带域 , 不再 随 时间而改变了 , 这就大大有利于数值计算的进行 (图 7 )。 这时 , 经过稍为仔细的换算 , 问题 (黔)就化成为如下的问题 : a戈 _ 1 0 么x 。 , , _ 一气二丁 ~ = ~ , - 气二尸气, 厂 , 二- ; 气 , U ‘‘、 U ‘、“ n O r (二竺一、一 o “- 、口哪 / 一 9 5 一 x l 。一。。 = 0 1 _ d x言 = 一日示 ,a : 1。一 。 x = S (t) ! “ _ 。 (28 ) S (0 ) = 0 这样方程虽然变成了非线性 , 但区域已为固定的了 。 文 献〔14 〕采用显式的差分方法求解问题 (28) : 戈梦+ ’一 x 梦 △t [型一护~ J = 1 , 2. 二 N 一 1 图 7 _ J乙么t 结果也很不错 。 他们还把此方法推广到二维问题 . 戈加 二 0 〔二:一喇飞概一 (2 9 ) 例如 , 对二维热传导方程 ut = u , + 为, (3 0 ) 如果把 y 看作是 二 , “ , t 的函数 , a g at 那么类似地可把 (30 》化成 : = _ 「旦丛. 一 一一 ,上一 旦兰芝1 一旦芝 LO x z ( 一旦丛、吕 。“ , J 。“ 、 a u l (3 1 ) 分界面边界条件 (7 )可以化成 : 斋 ·贵【, + (斋)’」· {“, (斋)一 ’ 一 k :(斋)一 ’} (3 2 ) 同样可以在(戈 , “ , t) 空间上一层层求解。 这个方法的好处是能把区域固定化 , 从而使求解大大简化 . 但有些问题 , 它不能计算 . 例如 : 问题 (27 )中的边界条件改为 : “】, _ 。 == f (t) = A 日in o t 这时 , 边界上温度开始时升高 , 后又下降 , 这样同一时间 , 不同的几点上可能出现同一个温度值 (即温度 “ 不是 (x , 川的单值函数 )时 , 方法即失效 。 2 . 贴体坐标系方法 它的思想是利用 自变量变换把运动区域变成为固定区域 . 例如如考虑上面的一维单相 S te fa n 问题 (2 7 ) , 二 的变化区域为 [0 , S (t )〕, 如果作如下自变量变换 : 七= 二 / S (t) (33 ) 于是 , 问题 (27) 变成了 : a u l 口Z u 乙 d ￡(t) _ _ 二二 _ _ - ~ - 二几- . ~一 .a t 5 2(t) a叠2 5 (t) d t“ = : 。 在 己= 0 S (0 ) = 0 贵 , ”< t<1 上 “ = 0 1 S (t) a“ _ _ q d s a七 一 尸 d t}在 息= 1 上 (34 ) 一 9 6 一 如果对 u 采用显式差分求解 , S (O , ds / dt 均可采用落后一次计算 , 那么就不必进行迭代。 S a it o h 〔‘. ]曾用类似的想法处理 了二维问题 。 例如考虑单相的 S te fa n 问题 。 (s a ) 、 (g b ) 、 (10 ) 、 (1 1 ) 、 (12 ) - 在固相域中心取一点 。 作为坐标原点 , 可 以建立一个极坐标系(r , 0 )。 设外边界 aD 的极坐标方程为 : r = r a (0) (3 5 ) 分界面的极坐标方程为 r = r r (0 , t ) (3 6 ) 设 耘二 C o n , t : 于是 , (s a )相 当于 琴 = 。v : T = 。 . 「李凛, ‘, 琴)+牛婴1口“ T 口 T \ 口 r , T 口 U Jr r ( r < r a下面只考虑方程(3 7 )而暂不涉及分界面条件 .如果引进如下的变换 : , _ r 一 介(0)七 = 一瑞‘杀黔生一不r r (0 , r ) 一 ra (9) 0 ( 0《2 究 (3 7 )r 一 r r (8 , 诊)r一(O) 一 r r (8 , t) (3 8 ) 那么 , 求解区域 r r ( r 《ra 0 《e( 2兀 就变成 了 : 0 ( 息( 1 (39 ) 0 ( 0夏 2北 这样 , 变区域就变成了一个固定区域(图 8 ) 。 这时方程 (3 7 )就变成 了一个相当复杂的非线性发展 _ _ J ~ 、 , _ ‘ , _ _ , 、 、 ,. 二、 . ⋯ _ 、 , ~ ~ , .⋯ _ . 、 . 一 、 , : _ , _ 。 一 ~ 一 、一 ~ 、、 , 、 ~ ~ a r , “ 、 , ~方程(见文献〔2 7〕) . 类似地 , 边界条件也可以变过去 . 如果采用显式有限差分法 , r , 和共子~的计算 /J ”工 、儿人私‘“ ’ 刁 ‘ . 入 ’叭 ~ ’~ , 「小 ” 。 ’ J M 人~ 升 . 州 / , 、月、 声 ’J ~ 少、 门 『‘、乙 , J ’“” 几 ” , at ” J 扩 ’一 落后一个时间步 , 那么计算就可以直接推进 . 文献〔1 9〕给出了不少有趣的结果 , 结果相当好 。 口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口口「「口口曰曰日日口口口口口口口口口门门口口 这个方法可用于多种问题 , 但对多连通区域的处 理 尚 有 困 难 , 显式时间步长较小等缺点 。 顺便提一下 , 贴体坐标方法最早 是由 A . M . W in s la w 【之? ] 提出的 , 它是数值网格生成方法的 一 种 , 目前已有较大的发展 [“吕] 。 3 . 引入热烩或热容量H 的办法【2“, ““, “. ] 〔20 〕曾引入了一个热容量H (H ea t C o nt e nt ) , 它是显热容量毛 ‘ “” 刁 目 J ’‘ “ J ‘ : ‘、‘ 一只 ~ “ 、 - - - - - - - 一 - - 一 ‘ ’ 目 ~ ~ “ 、‘ 曰 一 (S e n s ib le H e a t C e n te n t )和潜热 入之和 , 它与温度的关系式如图 9 所示 。 这样 , 方程式 (2 a) 和(Zb) 可统一成 : 李‘沐攀卜 p 擎 , 。< 二< ooU 荡 、 口苏 I 口不 (如 ) 在分界面 二 二 S (t )上 , 权 T )和H (T )有跳跃变化。 〔20 〕提出用一个有大梯度的光滑函数来逼近 无 和H , 这样分界面就 自动消失了 。 方程 <如)可化成如下非线性方程 : 一 97 一 a 八 _ , , 、一可一一 I 儿 气Z 夕 O 劣 、 旦Z止、日里 /aT aH= p we己了~ aTd t (41 ) 〔29 〕提出用折线来逼近图 9 的曲线 。 设 k 二 co ns t , 则H (T )可如下近似 : H (T ) = 几T , T 镇T , ‘ 。 H( T) = H( T一 。) +李 (少 一 T . + : ) , T , 一 。、T 、T . + 。汤e H (T ) = H (T 。 + 。) + C : (T 一 T一 e ) , 〔29〕采用完全隐式格式离散方程(40 ) : T ) T 。 + e (42 ) H 罗+ 且一 H 梦 p 一一一一万二 , 一一- 一 ” V ZT 梦+ ’ (43 ) 从而求解可比较稳定 。 这个方法有其一定的好处 , 不仅对一维是合适的 , 对高维也同样适用 。 但每一时刻必须解一 个非线性椭圆型方程的边值问题 。 可惜文献〔2幻解非线性方程组时只用到 G au : s一Se ide l迭代法 . 其它的 SO R 、 A D I方法等似乎更有希望 。 (三 )有限元方法 r , , ·“, · , 一l 由于分界面是移动的 , 因此通常的有限元法在 st e fa n 问题的处理上尚进行不多脚 , . L ync h & G r a y( 19 铭 , 19 80 , J . C o m p . Phys )曾发展了一种运动网格的有限元法 , 似乎是相当可行的 . 1 . 运动网格有限元法 由于边界 (单相问题 )是运动的 , 故节点自然也应是移动的 , 故应为时间 t 的函数 。 另外 , 有 限元插值函数小也应是 t 的隐函数 , 则有 : T (二 , t) == 艺币, (二 , t ) · 犷‘(t ) ~ 才 ~ (44 ) 子 一资鲁 ,仑 , , ) ·犷广票 如果对坐标 二 采用另一种插值函数 讥 , 则有 x = 艺 戈 ; (t) ~ 言 L y n e h & G r a了证明 T ; · 小‘ (x , t) (4 6 ) 一 98 一 (4 6 ) .由‘ V X”J 口 一一一票孔 a 于是有 : . , d 戈艺裂 」小, 一 = ,子 · v T“ 不 “不 (4 7 ) 这样 , 对热传导方程 : 。 OT , , _ 。 、 P七 ~花石, = V ’ 吸‘ V l 少 口‘ 可以用 G al e r k in 过程进行内积并分部积分 , 则可得 <小‘, C p 小, >争 一 <小‘ , p C 贵 · v小‘>T J < + 7 小‘, k V 小I ) T , 一 T , 丁。 n · (小r无V 小, ) d r = 0 (4 8 ) 其中 ( a) , ( b) 表示通常的内积积分 。 这样就得到了常微分方程组 (48 ) 。 对导数 d T , / dt 可用前差近似 d T , _ 刁了 ~ T 丁+ ’一 T 罗 t 。+ 一 t. 对 T , 项可采用显一隐格式 : T 一二 OT , 干 ’ + ( 1 一 0 ) T 二 ( 0 ( 0 ( 1 ) 从而 ( 4 8 )式可写成 : ( C奢, + OB梦, ) T 梦+ ‘ = [ C李, + ( 0 一 1 )B李,〕 · T 了 (4 9 ) 代数方程组 ( 49) 是非线性的 , 但若对 d 二 / dt 采用落后一次计算 , 方程组是线性的 。 他们的计算表明 , 采用线性插值函数不太好 , 而用 H e r m it e 插值则很好 [““〕。 2 . 变区域上的变分有限元法 E . V ar og hi 和 W . D . L . Fi n n 曾发展了一种变区域上的变分有限元法 , 成功地解决 了 一 维 单相 S t e fa n 问题 [“5 1 . 作者和 V a r o g lu 又把此方法加以推广 , 解决T 二维两相的 S t e fa n 问题t“‘l 。 该方法是先用线方法( M e th o d of Li ne )对时间 t 用后差进行隐式离散 。 化成每一时刻 的 椭 圆型方程的自由边值问题 , 从而化成等价的变区域上的变分问题 , 最后利用有限元法进行数值离 散 , 得到一组非线性方程组 , 从而可用 N e w to n 一R a p h s o n 方法进行迭代求解。 为了简单的说明该方法的过程 , 现以一维单相 St ef an 问题 (扬)式为例说明之 。 先 对 。, 项用 线方法进行隐式离散 : u 丹 一 u n 一 1 △t 一心二 O < 劣 < S (尹) (印 ) (为方便起见 , 以后均略去上标 “ 。 ” ) 引入所谓主泛函 ‘(一 , 二J:{晋·, +去 · (普一 ‘ )}d 二 一 I: F ‘一 , d X ( 61 a ) 注意 , 主泛函 I 不仅是未知函数 “ 的函数而且是区域边界函数 S 的函数 , 故是一个变区域上的泛 函数 。 对( 61 a) 式作一阶变分(它是变差的线性主部 ) , 一 9 9 一 各j = j(。 + 色。, 亏+ 乙: ) 一 I (。 , S ) 一 { ’“丫粤(“ + 。。) : +华 一 (。 + 。“) (共丝 一 “一 ; 、飞t ) 0 k Z 。 了 、 艺 , J dx 一买待, 十去 ·(普一 扩一 )卜 I {F ( : + 乙“ , u . + 各u 二 ) 一 F (。 , u 二 )}d 二 + {: ‘“ F ‘一, d ! 分部积分J: {- 、一去‘一 ‘,}“· + 二“·芡+ 【晋·, +去 “(昔一)}, “· 由于分界面上 故 “ 1 . = 0 f , r l , . _ 1 、 1 。 0 1 = 一 1 1 ““ 一 一飞二一 气“ 一 “ ‘ 少 1 0 “口 X J 0 L O r J 由于 各川。 一。 (, , 二 各以见图 10) . 即固定域内 “ 的变分 , 故由图上的几何关系知 : 0 < 义 < S (t ) 各“ 二 乙“ x = S (t)上 , 各“ = 乙。 + 。二各s 。. 色“ !。 。 , (, ) = “二 色n 二 “二 (各“ 一 u 二乙s ) = ‘加 一 “二阮 对对 , ‘砂沙汤尸 f ‘ f l , . _ , 、 1犷0 1 二 一 1 1 林“ 一 一气r 丫一 气“ 一 “ ’ 一少 1 0 材 口 X J o L 己了 J 1 , : 。 , 了一 百‘ }‘司 ”‘ 一叫 ‘一 。 。u0 ‘(t) 名+ d’ . 图 10 因此 : “‘+音S ‘“一 丽一买卜。 1 , , 一 百 诀‘一 “一月 ’OS 一 气如 十 l, ! ‘· 。 1 , _ 。 _ . 、 1 石二二万一 ’ 一 离丁- 气“ 一 “ 一1 10 ““ x乙二不 J 5“。 对照 (1动 式可知 , 如下变分为零的问题完全等价于 (1幻式的边值问题 . 。 v 一 。 7 . 1 2 泞” 一 ‘卜 ‘、: 。 犷 。 0 司八 ~ 0 1 十 、二 ~ I ~ es一- , r 二- - 一. 1 0 召 一 0 “n = U艺 、 凸犷 I ( 61 b ) 这就是变区域上的变分原理 。 对变分问题 (肚b) 的离散化 , 与通常的有限元相类似 , 但要注意的是目前不仅 “, 是未知数 , 夕 也是未知数 。 故采用通常的有限元后 , 可以得到 : I 二 I(“, , “ ,一“、 一 , , s , , △t , u f. 一 ’) 这里 “。⋯枷 一 : 为内点上的温度值。 因此 盯 = 乙 al 。 . al 一‘一 O 路 奋 十 一不 - 了 0 5 -d “j ‘ d 了’ 一 100 一 从而可得方 + 1 个非线性代数方程组 : f , (u 。 , “: ⋯ “, 一 , , S ” , △t , : 卜 ‘) 二 o (j = o ⋯N ) 我们可用 N e w to n 一R a p h s o n 方法迭代出 产 时刻的温度分布 : 佘二 u:- _ ; 及运动边界的位置 S 气 这个方法的好处在于把所有未知函数 (温度和运动边界位置 )同等看待 , 把所有边界 条 件(包 括分界面上两个条件 )同时包含在变分问题之中 , 因此精度较高 。 这种做法相当于有限元法 中 的 混合方法 。 这种方法的优点 , 除了精度高外 , 对各种边界条件均可应用 (包括第一 、 二 、 三类及辐 射 边 界条件 ) , 对维数似乎是没有限制的文献 (〔触〕解决了二维两相 5 te fan 问题 ) 。 它的缺点是 : 时间 离散后的方程必须存在相应的泛函 , 若不存在这种泛函就不能应用 . 另外 , 计算时间似乎较长 . 四 、 更复杂的 Stef a n 问题及边界元方法应用的设想 本节介绍某些更复 杂 形 式 的 S tef a n 问题 , 同时对边界元方法的应用谈一些初浅的 想 法 。 1 . 对流扩散 S t efa n 问题 前面各节所谈的相变传热问题 , 都是假设固相和浓相的密度 p . 和 p : 是 相同的 , 即 p : = p , = p 。 但是 , 实际情况并非如此 , 即使对冰和水密度也不一样 , p 冰/ p 水 = 0 . 9 7 。 当然 由于差别较少 , 一般可不予考虑。 但是 , 当密度差较大时 , 往往由于体积的变化会在液相区域产生一 种 对 流 现 象 , 这时液相区域的运动方程不再是普通的热传导方程 , 而是一个对流扩散方程 . 例如一维杆的固结问题 , 在 t < 0 时杆(0 < x < co )处于均匀 温 度 T ,( > T动 。 在 t 二 0 时 , 在 二 = O面上突然降温至常温 T 试 < T衬 . 这样一个固结面 二 = S (t )就向右移动了 (图 11) 。 设 P. > p ; 则有 : T . < T . p . e .擎 一李(、:擎、 。< 二 < s (, )O 犷 口 X \ O X / _ 。 了aT , . _ _ aT , 、_ a I t _ OT , 、 p ‘灿云火.万不 宁 “ ‘万王一) 一 一万犷灭凡 , 泊百-) S (t) < x < co 初始条件为 : 5 (0 ) = 0 T , (劣 , 0 ) = T , 0 ( 另 < co (石2 ) 边界条件为 : T 。 (0 , t ) = T . T x (S )r) , r ) = T . (S (t) , aT , . aT , 瓦. 竿三一 一 k i竺李上 = 入p :一 a x ’一 ‘ a x ’ , ~ t) 二 T . d s (t) J t 其中对流速度 晰 一 (色汁)票 (碍 ) 一般说来 , 由于 p : 一 p , 比较小 , 这种小对流速度问题数值计算没什么困难 . 但是 , 〔1〕中曾举出 如陨石的表面气化问题 , 这时控制方程为 : 一 10 1 一 Or . a T . ar . ar _ , , ~万了 十 “‘ ,刁王一 十 “f 百万冲 “ · ,百矛 一 “V 一 ‘ (以) 当对流速度 }厂 !相当大时 (即 Pec le t 数 尸。》 1) , 许多数值计算会产生较大的振荡和误差 , 必须要 用特殊的办法才行 t“’】。 这类间题在对流一扩散的情况下同样存在 。 2 . 合金的熔解 前面讨论的熔解 (或固结)问题 , 都认为在一个熔点温度下材料由固相转变成液相 。 一般说来 对完全纯的物质是这样的 。 但是 , 对于像合金这类材料 , 由于熔点不是真正的一个点 , 而是一个 熔解范围 (或区域 ) , 材料在这个温度范围内才逐渐的熔解 。 这种问题比经典的 St ef an 问题 更 为 复杂 . C ro w le y & O c k e n d o n 曾对合金的固结问题进行了数值研究〔““] , 他们采用的是热容 量 或 烩的方法 , 结果相当好 . 其它工作似乎不多 。 在应用上这类问题是值得引起重视的 。 3 . 关于应用边界元方法的问题 应用边界积分方法来解决 St e fa n 问题似乎已有较长的一段历史了 [ ’, 之] 。 但是 , 应 用 边 界元 方法来数值求解 S te fa n 问题的文章尚未多见 , 特别是高维问题 。 例如 , 考虑半无限液体的固结问题 (6 2) , 令 p . = p : = p . (b : = c : = c) , 那么 , 利用位于 工 = S (t) 上的运动热源的概念 , 题汇, , 2 2 , 、几 , 。 , _ , _ , _ _ 七 ,汉花乙 1 佑 = U , 凡 1 = 凡2 = ‘ , “ =—PC问题 (6 2 )可以化为如 下的 等价 I’d a Z T a x Z + 一 l d s 。 , 。 , . 、 、 1 a T , , _ _ _ p 八 一万: - . 0 气X 一 O 气r 少) 二一 一丈二一 U 《气X ‘、 ,即a 了 a o 了 T (0 , t) = T . = 0 T (戈 , O) 二 T , lim T (x , t) = T ‘ (56 ) T (S (t) , t) = T . 其中 各(二 )为 D ira c d el ta 函数 。 容易证明 (对方程式两边关于 二 从 S (t) 一 。 积分至 S (t ) + 。 , 然 后会 。。0) 〔Z J, 分界面上另一个条件 (6 )会自动满足的 . 为了求解 (5动 , 可令 T (工 , t ) = T i (x , t ) + T Z(戈 , t ) 其中 T , (x , O 满足如下方程和初边值问题 : a Z T I a劣 2 1 aT , = — 一一丈二一 X 碑) U a 口了 T (o , t) = 0 T (x , o ) = T l t> 0 (66 ) 戈 > 0 而 T : (工 , t) 满足 a Z T : , 入p d s 、 , ‘. c , 二 、、 - 万妥丁 , 一反厂 一万了U “涌 一 J 、‘ 户产 - 1 aT , _ 。 — .一砚二 X 剖夕 Ua O t T (0 , t ) 二 0 T (x , 0 ) = 0 t> 0 x > 0 (盯 ) 一 10 2 一 显然 , 问题 (石6 )的解为 t’, 艺〕 : :。二 , , , 一 : ‘e·f (厄箫) (印 ) 而 T : (x , t) 可用 G r e e n 函数来表示 , , . 、 a f云 f“ 。 , , . 、 , d ‘(T ) 。 , , , 、 , , 1 2 气劣 , 犷) = 二一 ! l 行 气x , r ; x ’ , T 少 · p 八一一不一 0 气戈 一 S 吸T 少J a x u TK J O J o U 不 其中‘(二 , t ; 二 ’, 口是半无限区域热传导方程的 G re e n 函数 (齐次边界条件) : ‘(x , 1; 戈 ’ , T ) 二 二一牛开 于气万一一书一艺、 / 究仅 吸r 一 T 夕 (x 一 义 , )” 4 a (t 一 下 ) 一 e X p ( 戈 + 戈 , )“ 4 a (t一 下 ))」 (6 9 ) 关于 x ’积分后可得 : ~ , 、 久 f, d s (T ) l r 1 2 气工 , 才声 = 不石一了节资产 ! 一〕了一, . 一万 二于 ! eX P艺七 , 犷 卫“ J O 二 V ‘ 一 T L (x 一 s (T ))2 4 a (t 一 T ) Td 一...‘J、,.产 一 e x P (戈 + S (r )) 2 4 a (t一 下 ) (6 0 ) 最后 , 由分界面条件 T , = T : (S (t ) , 才) + T : (S (t) , 矛) 可以得到 S (t) 的积分微分方程式 : , 。 。 了 S (t ) \ _ 入 (, d s (T ) 1 T . 二 T ‘ . e rf . 止二乏里注 】+ 二宁分之.六石 l 二鉴召一 尸一立一-‘ , 一 ‘ ’ 一 ““ 、2必丽j 「 ZC , (二a )’IZJ。 d t 杯才二万 「_ _ _ l (S (t) 一 s (T ) 2 、 _ _ _ 了 (S(t) + s (T )):、1 _ , _ _ , 。 一 、 1e x p l 一竺毕 ~共下竺共二一 e x p 一 望岑 ,钾二二 兰之己z 一 11 d t (6 1 )L 一 ‘ 、 4 a (t 一 丫 ) , 一否 、 4 a (t 一 T ) 月 一 ’一 ’ 对一维问题 , 可以利用S (t) ~ 记下的结果来得到比例系数的超越方程 , 这与精确解相一致 . 在作者所阅的文献内似乎尚没有人应用边界元 (线性 、 二次元等)的思想对 (6 1) 进行过详细的 讨论和计算 , 这似乎是很遗憾的 。 对高维问题 , 类似地可以应用运动热源的想法 , 把分界面 两边的方程统一成等价 的 一 个 方 程 : , “ T +粤。, 。(、 一 , 。) = 粤琴 在 D 二 + 刀‘内凡 , 认 口奋 (62 ) 其中 v 是任一点尸至分界面的法向距离 , v 。为分界面上点的位置 . 例如 , 当分界面方程为 z 二 S (戈 , 夕, r) (6 3 ) 时 , 方程可化为 : 护T 十粤 丝二兰竺‘典业。(二 一 : 。) 一李擎凡2 0 不 “ O 丁 (分连) 类似地可得到分界面位置的积分一微分方程式 。 由于边界元方法能使问题的维数降低一维 , 因此是一个很有吸 引力的数值方法 , 可惜目前做 的人很少。 当然 , 由于边界元方法得到的是满阵 , 求解不很快 , 这是它的最大缺陷 , 希望搞计算 传热的同志对此有所注意 . 一 加3 一 参 考 文 献 t l〕 e a r sla w , H . s , a n d Ja e g e r , J . e . , e o n d u e t io n o f H e a t i n s o lid s . : e 。 o n d e d . , Cra r e n d o n Pr e s , , L o 且d o 几 (195 9 ) . 〔 2 」 O z is ik , M · N · , H e a t C o n d u e t io n , J o h n W ilo y a n d S o n s , N e w Y o r k (19 80 ) . 〔3 〕 C a n n o n , J · R . , “M o v in g B o u n d a r y P r o ble m s , ” D · G · W ils o n e t a l; e d · , A e a d e m ie Pre s s . Io n , d o n . PP . 3一2 4 (1 9丫8 ) · 〔4 〕M o y o r , G · H · , 5 1人M , J· N 住 m e r · A n a l· 10 , 52 2一 53 8 (19 73 ) · 〔5 〕 O r n