在自然界和工业生产中，广泛存在着带有相变的热传导问题，例如：地球两极冰块的熔化和凝固、晶体的生长、铸件的固化等问题。这类问题的特点是区域内存在一个随时间变化的两相界面，在其两相界面附近发生相变，释放或吸收大量热量。相变材料这一特征可被用于储存能量或达到控制环境温度目的。
由于相变传热问题存在一个移动边界，因此相变热传导问题属于非线性问题。现有的求解相变热传导问题的数值方法主要有：有限元法、有限体积法、无网格法、边界元法等方法；对相变界面位置的处理方法主要分为两类：固定区域法和界面追踪法。如：Chessa等[1]结合扩展有限元法和水平集方法对凝固问题进行了研究；Zhang等[2]采用有限点法和焓法相结合模拟了金属凝固过程；Muhieddine等[3]采用有限体积法对比了潜热源法和显示热熔法处理一维相变问题的差异；李海梅等[4]采用等效热熔法推导了平面相变传热问题的有限元解；Yang和He[5]利用基于Sigmoid函数光滑化的等效热熔法和无网格伽辽金法建立了求解相变传热问题的数值模型。固定区域法不需要连续地追踪移动界面位置和网格重构，整个计算区域建立统一的方程求解，但是在处理等温相变问题时，必须人为假设一个虚拟的相变温度区间，相变界面位置必须采用焓或者温度的线性插值得到。而界面追踪法连续追踪相变界面位置，具有较高精度，无需插值计算。同时已有的相变热传导问题研究在时域上大多采用差分法处理，计算结果对时间步长较为敏感。
边界元也是常用的数值方法之一，如：王海涛和姚振汉[6]将快速多极边界元法应用到大规模传热分析中；Liu和 Lu[7]用双重互易边界元法求解了一维半无限大平板非等温相变问题；2002年Gao[8－9]提出径向积分边界元法，并用于求解瞬态热传导问题。Jesus等[10]采用径向积分边界元法分析了各向异性材料壳的动力学问题；邓琴等[11]采用径向积分边界元法分析了弹塑性静力学问题；Yu和Yao等[12－13]将径向积分边界元法与精细时域展开算法相结合求解了非稳态傅里叶热传导问题与非傅里叶热传导问题。边界元法具有求解的温度和热流密度同阶精度的优势，同时只需要在边界上划分单元，而界面追踪法能精确计算相变界面的位置。边界元法与界面追踪法相结合具有高精度、计算量小和能直接获取界面位置，无需插值的优势。
本文利用径向积分边界元法给出了单相凝固问题的边界积分方程，并结合精细积分算法[14]和界面追踪法[15]建立了单相凝固问题的数值模型，并通过两个算例，对本文方法进行了验证，得到了满意的结果。
1 基本方程
讨论如图1所示的平面单相凝固问题，其区域设为Ω，固定边界为O B∂，初始时刻域内含有温度为熔点温度m T的液体。在0t＞时，对域边界O B∂施加一低于熔点的温度或有热流输出，液体开始向内凝固，其中相变界面记为I()B t∂。假设自始至终液相区温度保持不变，则问题为单相凝固问题，只需要分析固相域Ω即可。
直角坐标系 1 2Ox x下固相域问题的控制微分方程为：
式中：k为热导率；b为热源项；ρ为材料密度；p c为比热容。注意：本文采用爱因斯坦求和约定，即相同下标表示求和。
由于k为常数，可引入规格化温度
和常数
则控制方程可改写为：
相对应问题的初始条件为：
而固定边界O B∂上的边界条件可假设为：
式中：
图1 凝固问题的示意图
Fig.1 Geometry of solidification problem
在相变界面I()B t∂上始终有：
同时因为相变问题会在相变界面I()B t∂上产生热交换，其能量平衡方程可以写为：
式中：L为相变材料的凝固潜热； n V为相变界面移动的法向速度。
2 时域径向积分边界元法
对控制式(2)本文采用边界元法求解。对 mt时刻，首先引入权函数u*，给出控制式(2)的加权余量式：
对式(8)左端第一项域积分进行两次分部积分，并采用高斯散度定理有：
式中，边界
取u*为平面稳态热传导问题的基本解，即有：
式中，r为源点 p x和场点 q x之间的距离，于是：
这里当 p x位于域内部时 1c=；而当 p x位于边界时
α为源点两侧边界的夹角。
将式(11)代入式(9)，可得：
再将式(12)代入式(8)，整理有：
式中，
式(13)中仍含有域积分项，本文采用径向积分法[8]将其转换为边界积分。
对于式(13)中已知函数的域积分项
可以直接采用径向积分法将其转换为边界积分。
式中：
x Q 表示边界上的场点； x q表示计算域内的场点，可以是内部点，也可以是边界点；ξ (x p , x q )表示源点 x p和场点 x q的距离，如图2所示。
图2 点 p x、 q x、 Q x和积分变量ξ的关系
Fig.2 Relation between points p x,q x,Q x and the integral variable ξ
而对于未知函数/T t∂∂~，不能直接使用径向积分法。首先采用四阶样条紧支径向基函数和附加的线性多项式对其近似：
式中，
为径向基函数：
N为节点总数，包括所有的边界节点和内部配点
表示场点 q x与当前节点坐标 i x之间的距离， id为当前节点的紧支半径。附加的线性多项式中存在3个未知量，需要添加3个补充方程[16]：
它和全部N个节点处温度的时间导数相等条件联立可得到方程组：
其中：
当任意两节点坐标不重合时，Φ为非奇异矩阵，则有：
将式(16)和式(21)代入式(13)的最后一项域积分中，采用径向积分法[8]整理得：
式中：
于是式(13)转换为一个纯边界积分方程：
假设边界∂Ω离散为 N E 个线性单元，共包含N b 个边界节点，在给定温度边界条件Γ1上分布N b1个节点，而给定热流密度的边界条件Γ2上分布N b 2个节点，且
同时，在域内有 N I个配点。注意本文下标b代表边界量，I代表域内量。则待求的未知量有给定温度边界上的节点热流密度 q b 1、给定热流密度边界上的节点温度 T b 2和域内配点温度 T I，令
将式(24)中的场点x q分别取为边界Γ1上的N b1个节点、边界Γ2上的 N b 2个节点、和域内 N I个配点，则有：
式中： b1A 、 b1A 和 b1A 为已知量，是由已知变量
和热源项得出。
对式(25)中的系数矩阵整理有：
由式(26)第一式可以得：
将式(27)代入式(26)第二式整理可得：
其中：
根据常微分方程组的求解理论，式(28)对应的通解可以表示为：
式中，
为矩阵指数函数。
求解式(30)即可给出 mt时刻瞬态问题(2)的数值解，并给出相变界面
上的热流密度。再根据式(7)计算出 mt时刻节点移动的速度向量 m=V
其中M是相变界面上的节点总数。
式(30)可采用精细积分算法[14]求解，该方法是一种显示积分法并且是无条件稳定的。这里需要说明的是，式(30)最后一个积分项中的被积函数F是由已知的边界条件和热源项生成的。当问题边界条件和热源函数为多项式、指数函数和三角函数等函数时，该项积分可以解析求得。而对于一些较复杂的情况，可在时段上采用线性插值获得其数值结果。
3 界面追踪法
依据mt时刻相变界面I()B t∂上的热流密度，就可以采用Zabaras和Mukherjee[15]的界面追踪法给出相变界面移动的位置。
图3所示，利用包含界面节点i的两个相邻单元的单位外法线向量 i n和 1i-n 来确定节点i移动的方向：
式中， il和 1il-是第i和第 1i- 个单元的长度。
显然方向 ()i n 的单位向量为：
同时，利用获得的mt时刻相变界面I()B t∂上的热流密度，通过式(7)确定变时间步长1mt+Δ：
图3 节点移动的方向
Fig3 Direction of node moving
式中，sΔ是每个时间段给定的最大凝固厚度。
然后是确定相变界面节点 1mt+时刻的移动速度及相变界面节点移动的位置。因为当前问题属于非线性问题，为保证求解的精度，需要采用迭代法[15]进行求解，其
时间段的主要步骤如下。
1) 根据 mt时刻相变界面上节点的移动速度
先给出 1mt+时刻一个预测节点移动速度向量
2) 利用当前 1mt+时刻节点移动速度向量
就可给出相变界面上节点的移动位置：
其中，
代表在mt时刻动边界上第i个节点的位置。即获得 1mt+时刻新的动边界
3) 采用第 2节的时域径向积分边界元法分析
时刻瞬态方程式(2)，获得新的 1mt+时刻的相变界面上节点的移动速度 V~ m +1。
4) 收敛性检测。如果条件
满足，即认为
时段的计算已经收敛。否则，用
返回第2)步。
需要注意的是，在更新相变界面的节点位置后，要检查相邻节点之间的距离。对于向内凝固问题，相变界面会越来越小，可能出现相邻节点距离太近，此时产生近似奇异性，造成计算结果不稳定甚至失效，需要对相邻太近的节点做移除操作。而对于向外凝固问题，相变界面会越来越大，需要适当的增加单元节点，同时为了保证计算精度，当凝固区域达到一定的大小时，需要增加一定量的内点。
4 数值算例
算例 1：分析图 4所示的变温圆管。设圆管半径
圆管表面温度为给定温度：
式中：指数积分函数
圆管内介质的热传导率
初始时刻整个区域处于未冻结状态，其温度等于相变温度即
图4 算例1初始时刻的求解模型
Fig.4 Computational model at initial time for example 1
该问题属于单相向外凝固问题，所以只需考虑冻结区域。本问题径向冻结位置的解析解[17]为：
为了程序便于启动，假设圆管已经冻结了0.1 cm的冻结层，此时
取最大凝固厚度1 cm sΔ=，收敛误差
规定凝固厚度每增加10 cm，增加20个内点。表1给出了采用本文方法所获得的不同时刻移动界面径向位置，其与解析解非常吻合。求解过程中，每个时间段的迭代次数均不超过 5次，当
时，迭代次数均不超过3次，收敛速度很快。
表1 移动界面径向位置的误差
Table 1 Errors of locations of freezing interfaces
时间/(×105 s) 本文解/cm 解析解/cm 误差/(%)0.475 003 4 6.746 480 79 6.713 066 65 0.497 747 74 1.615 965 8 12.457 933 62 12.381 938 60 0.613 757 07 3.475 442 0 18.272 067 18 18.158 400 21 0.625 974 59 6.063 115 5 24.136 572 48 23.983 937 52 0.636 404 97 9.382 671 7 30.029 406 49 29.835 671 96 0.649 338 58 14.184 411 3 36.927 826 95 36.684 129 40 0.664 313 32 19.098 961 2 42.854 500 83 42.567 431 09 0.674 388 22 24.751 847 4 48.790 279 65 48.459 214 55 0.683 182 96
算例2：如图5所示，考虑一个2m 2m×的正方形区域且不带热源的凝固问题，其中热传导率
密度
比热容
潜热
初始时刻整个区域处于未冻结状态，其温度等于相变温度即
而方形表面为给定温度边界条件：
图5 算例2初始时刻的求解模型
Fig.5 Computational model at initial time for example 2
该问题是一个单相向内凝固问题，即只需考虑固相域。本算例在计算前期远离角域的相变边界是平行于坐标轴的，所以文献[15]的半无限角域的解析解可以作为本文的一个很好的近似解。给定一个初始的凝固厚度为0.1 m的区域，此时 0t=0.00939 s。
取最大凝固厚度
规定x轴上凝固厚度增加0.2 m，域内增加40个内点，同时当相变边界上单元长度变为原来长度的1/5时，将单元的两个节点合并到单元中间位置成为一个节点。表 2给出了采用本文方法所获得的不同时刻移动界面在x轴上的位置，其结果与解析解非常吻合。图6给出了同一种单元划分情况下，不同最大凝固厚度Δs (即不同时间步长)时x轴上相变边界位置的相对误差。从图中可以看出Δs分别取0.03 m、0.04 m、0.05 m和0.06 m时计算结果仍然有很好的精度，结果对时间步长不敏感。
表2 凝固界面在x轴上位置
Table 2 Locations of freezing interfaces on x-axis
时间/s 本文解/cm 解析解[15]/cm 误差/(%)0.012 978 24 0.118 997 09 0.117 768 10 1.043 563 52 0.034 999 99 0.192 783 62 0.193 398 76 0.318 069 21 0.064 580 68 0.261 445 76 0.262 706 62 0.479 949 48 0.093 727 83 0.318 128 51 0.316 485 63 0.519 100 46 0.127 766 12 0.370 758 12 0.369 511 18 0.337 457 67 0.172 714 98 0.430 715 45 0.429 620 23 0.254 925 61 0.209 095 68 0.474 436 64 0.472 707 24 0.365 849 71 0.277 143 65 0.549 404 57 0.544 217 11 0.953 198 01
5 结论
本文针对平面单相凝固热传导问题，提出了时域精细积分边界元方法与界面追踪法相结合的数值
图6 对于不同的sΔ, x轴上凝固界面位置的误差Fig.6 Errors of solidification front on x-axis for different sΔ模拟方法。数值算例结果表明，本文方法具有迭代次数少、求解精度高、计算结果对时间步长不敏感的优点，是求解相变热传导问题的一种有效数值方法。本文扩展了径向积分边界元法的应用领域，同时也为径向积分边界元法求解复杂的相变问题奠定了基础。
[1] Chessa J, Smolinski P, Belytschko T. The extended finite element method (XFEM) for solidification problems [J].International Journal for Numerical Methods in Engineering, 2002, 53(8): 1959－1977.
[2] Zhang L, Rong Y M, Shen H F, et al. Solidification modeling in continuous casting by finite point method[J]. Journal of Materials Processing Technology, 2007,192(5): 511－517.
[3] Muhieddine M, Canot E, March R. Various approaches for solving problems in heat conduction w ith phase change [J]. International Journal on Finite Volumes,2009, 6(1): 1－20.
[4] 李海梅, 顾元宪, 申长雨. 平面相变热传导问题等效热容法的有限元解[J]. 大连理工大学学报, 2000,40(1): 45－48.Li Haimei, Gu Yuanxian, Shen Changyu. Finite element solution of heat transfer w ith planar phase change by equivalent heat capacity method [J]. Journal of Dalian University of Technology, 2000, 40(1): 45－48. (in Chinese)
[5] Yang H T, He Y Q. Solving heat transfer problems w ith phase change via smoothed effective heat capacity and element-free Galerkin methods [J]. International Communications in Heat and Mass Transfer, 2010, 37(4):385－392.
[6] 王海涛, 姚振汉. 快速多极边界元法在大规模传热分析中的应用[J]. 工程力学, 2008, 25(9): 23－27.Wang H T, Yao Z H. Application of fast multipole boundary element method on large scale thermal analysis[J]. Engineering Mechanics, 2008, 25(9): 23－27. (in Chinese)
[7] Liu J, Lu W Q. Numerical simulation of non-isothermal phase change problem using a DRBEM w ith augmented items [J]. Journal of Engineering Thermophysics, 2006,36(7): 408－416.
[8] Gao X W. The radial integration method for evaluation of domain integrals w ith boundary-only discretization [J].Engineering Analysis w ith Boundary Elements, 2002,26(10): 905－916.
[9] Yang K, Gao X W. Radial integration BEM for transient heat conduction problems [J]. Engineering Analysis w ith Boundary Elements, 2010, 34(6): 557－563.
[10] Jesus L J M, Cimini C A, Albuquerque E L. Application of the radial integration method into dynamic formulation of anisotropic shallow shells using boundary element method [J]. Key Engineering Materials, 2015,627(9): 465－468.
[11] 邓琴, 李春光, 王水林, 等. 无域积分的弹塑性边界元法的非线性互补方法[J]. 工程力学, 2012, 29(7): 49－55.Deng Qin, Li Chunguang, Wang Shuilin, et al. Nonlinear complementarity approach for elastoplastic bem w ithout internal cell [J]. Engineering Mechanics, 2012, 29(7):49－55. (in Chinese)
[12] Yu B, Yao W A. A precise time-domain expanding boundary-element method for solving three-dimensional transient heat conduction problems w ith variable thermal conductivity [J]. Numerical Heat Transfer Part B Fundamentals, 2014, 66(5): 422－445.
[13] Yu B, Yao W A, Zhou H L, et al. Precise time-domain expanding BEM for solving non-fourier heat conduction problems [J]. Numerical Heat Transfer Part B Fundamentals, 2015, 68(6): 511－532.
[14] 钟万勰. 子域精细积分及偏微分方程数值解[J]. 计算结构力学及其应用, 1995, 12(3): 253－260.Zhong Wanxie. Subdomain Precise integration method and numerical solution of partial differential equations[J]. Chinese Journal of Computational Mechanics, 1995,12(3): 253－260. (in Chinese)
[15] Zabaras N, Mukherjee S. An analysis of solidification problems by the boundary element method [J].International Journal for Numerical Methods in Engineering, 1987, 24(10): 1879－1900.
[16] Gao X W. A meshless BEM for isotropic heat conduction problems w ith heat generation and spatially varying conductivity [J]. International Journal for Numerical Methods in Engineering, 2006, 66(9): 1411－1431.
[17] O'Neill K. Boundary integral equation solution of moving boundary phase change problems [J].International Journal for Numerical Methods in Engineering, 1983, 19(12): 1825－1850.