D0I:10.13374/i.is8n1001053x.1991.03.023 北京科技大学学报 第13卷第3期 Vo1.13No.3 1991年5月 Journal of University of Science and Technology Beijing May 1991 一种采用热焓处理相变导热 问题的有效算法 何评·俞昌铭· 摘要：通过对采用热培处理相变导热向题的有限元分桥，提出了消除选代求解的一种有 效算法。该法克服了传统算法的缺点。从计算结果看，该算法具有精度高、机时少的特点。 关键词：培，有限元，数值，导热，相变 An Effective Numerical Method of Heat Conduction with Phase-Change by Using Enthalpy He Ping'Yu Changming' ABSTRACT:According to the finite element analysis of heat conduction with phase-change by using the enthalpy,an effective numerical method to eliminate the solution of iteration is presented.It overcomes the disadvantages of the former numerical methods,and it is shown from the new numerical results that this method has the characteristics of high accuracy and needs less calcu- lating time。 KEY WORDS:enthalpy,finite element,numerical,conduction,phase change 处理相变导热问题的方法有许多，其中，焓法是数值求解这类问题最常用的方法之一。 由于这种方法把格和温度一起作为待求函数，控制微分方程中包含了两个待求变量，所 1990一11一12收稿 ·热能工程系(Department of Ene rgy Engineering) 228
第 卷 第 期 ， 年 月 北 京 科 技 大 学 学 报 盯 。 ， 一种采 用热烙处理相变导热 问题的有效算法 何 评 ’ 俞昌铭 ， 摘 要 通过 对采 用热始 处理 相变导热问题的有限元分析 ， 提 出 了消除迭 代求 解的 一种有 效算法 。 该法克服了传统算法的 缺点 。 从 计算结果看 ， 该算法具有精度高 、 机时少的特点 ， 关健词 烤 ， 有限元 ， 数值 ， 导热 ， 相变 一 夕 ” 夕阴 刀 夕 一 ， ， 几 ， ， ， ， 处理相变导热 问题的方法 有许多 ， 其 中 ， 烤法是数值求解这类 问题最常用 的方法之一 。 由于 这种方法把烙和温度一起作 为待求 函数 ， ，石控制微分方程 中包含 了两个 待 求 变 量 ， 所 一 一 收稿 热能工程 系 DOI ：10．13374／j ．issn1001－053x．1991．03．023
以，离散后最终得到的n个代数方程包含着2个变量，从而使代数方程组的求解必须采用 高斯一塞德尔迭代或与其相似的异步迭代：1。这种异步迭代收敛速度很慢，尤其与相变的 非线性特性交织在一起，每个时间步通常需要很多次迭代，很费机时，因而，计算起来很 不方便。为了解决由于迭代所带来的缺陷，本文用有限元对相变导热的焓法进行了研究，提 出了一种有效的算法，即通过对单元热导矩阵的巧妙处理，可将非线性导热中采用的三时间 格式引入焓法方程的求解，从而不需要送代就能得到解的结果。因此，可大大减少计算工作 量，使计算时间大为减少。 1 有限元半离散方程 一般说来，相变导热培法的控制微分方程通常可写为： (pH) ot =div(kgradT) (1) 其有限元半离散方程为： 〔S〕rP)+〔K〔T)=〔F) (2) s-。pM,d2 (3) k,=∫MvM,d0 (4) 1.-f qMar (5) 上式中，Q是计算区域，P是计算区域的边界，M:,M;是节点位置，〔S门是密度矩阵 〔P〕是姑矩阵，〔K)是热导矩阵，〔F门是热源矩阵，〔9〕是源项。从上述有限元半离散方程 可以看出，方程(2)还不能直接求解，需要进一步离散。 2 时间离散 如果引人三时间格式离散不稳定项，则方程式(2)可变为2： CS)〔Pa】"+1+〔K){(0.5+a)〔T)"+1÷(0.5-2a)〔T) +a〔T)-1}=〔F) (6) 其中，a是任取常数，在0~1范围内变化，另外： CPH"+1=(CH)+1-〔H))/△t (7) 〔T)+1=〔T]+△〔P)+1 (8) 上标*表示矩阵和向量取t时刻的值，即： 229
以 ， 离散后最终得到的 打 个代 数方程 包含着 个 变量 ， 从而使代数方程组 的求解必须采 用 高斯一塞 德尔迭代或与其相似的 异步迭代 〔 ‘ ’ 。 这种异步迭代收敛速度很慢 ， 尤其与相 变 的 非线性特性 交织在一起 ， 每个时间步通常需要很多次迭代 ， 很 费机时 ， 因而 ， 计 算 起 来 很 不方便 。 为了解决 由于迭 代所带来的 缺陷 ， 本文用 有限元 对相变导热 的治法进行了研究 ， 提 出了一种有效的算法 ， 即通过对单元热 导矩阵的 巧妙处理 ， 可将非线性导热 中采用的三时间 格式引入活法方程的求解 ， 从而不需要迭代就能得到解的结果 。 因此 ， 可大大减少计算工作 量 ， 使 计算时 间大 为减少 。 有限元半离散方程 一般说来 ， 相变导热焙法的控制微分方程通常可写 为 舟 其有限元半离散方程 为 〔 〕 〔 〕 〔 〕 〔 〕 二 〔 〕 ， 二 歹 。 ， ‘ 了‘ “ “ ‘ ， 于 。 ‘ ‘ · ， “ ， · ， 。 “ 上式 中 ， 口 是计算区 域 ， 是计算区域 的边界 ， ‘ ， ， 是节点 位置 ， 〔 〕 是密度 矩 阵 〔 〕是治矩 阵 ， 〔 〕是热导矩阵 ， 〔 〕是热 源炬阵 ， 〔的 是源项 。 从上述 有限元半离散 方 程 可 以看出 ， 方程 还不 能直接求解 ， 需要进一步离散 。 时 间 离 散 如果引人三 时间格式离散不 稳定项 ， 则方程式 可 变为 〔 ’ 〔 〕 ， 〔 二〕 ’ ‘ 一 ’ 〔 〕 。 〔 〕 ’ 十 ’ 一 〔 〕 ’ 〔 〕 一 ‘ 二 〔 〕 ， 其 中 ， “ 是任取常数 ， 在 。 范 围内变化 ， 另外 〔 二 〕 ” 十 〔 〕 ” ’ 一 〔 〕 “ 〔 〕 ” 〔 〕 ’ △考〔 〕 ’ 十 上标 带 表示矩 阵和向量取 时刻 的值 ， 即
”=t.+女或”=. 对应t有： T=1.5T"-0.5T-1或T*=T (9) 若以格的变化率为待求变量并将式(8)代入式(6)则可得： CS〔PJ+1+CK)(0.5+a)△〔P)+1=〔F)'-〔K){(1-a) 〔T)"+a〔T)-1} (10) 对于格和温度而言，它们的变化率存在如下关系： 若无相变发生，则 〔PH)=〔C〕〔P) (11) 因而有： 〔PH+1=〔C]+1〔PJ+1 (12) 其中，〔C)和〔C)+1是比热矩阵。 若有相变发生，对于单元物质的相变向题来说（即无二相区），由于正在发生相变的节 点其温度不发生变化，故(11)、(12)两式不再适用。但若采用下述方法处理相变节点所在单 元的热导阵〔K),仍可使两式成立，例如，如果某单元的第个节点处于相变，则令该单元 〔K)中第行和第列中的元素为零，：=0意味着可以维持节点的温度不变，因此，将 关系式(11)代入方程(10)可得： {CS)·+(0.5+a)AtCK'/CC)}CPa)"*1 =〔F)'-〔K〕·{(1-a)〔T)+a〔T]-1} (13) 利用方程式(9)，（本文在实际计算中采用T·=1.5T"-0.5T-1的方法）； 求解方程式(13)可得〔P门+1，然后用式(7)求出〔H)+1,再由培和温度的关系曲线可以求 出CT)+1。 由上面的推导可以看出，引入三时间格式以后，可以不需要迭代而达到计算一次完成的 效果。 了数值计算 为了说明上述数值计算方法的可行性，本文就一维半无限空间的凝固问题和二维无限大 方形角区液体的凝固问题进行了计算，在计算中，采用双线性四边形等参单元，单元热导阵 的计算采用补偿后的一点高斯积分〔3？4)，同时，还采用了集中质量法处理质量矩阵的方 法。 图1所示是本文采用的一维半无限空间的简化模型和有限元网格划分图。为了便于比 较，本文采用如下的物性数据，即： 230
‘ · 十 或 ’ “ · 对应 有 ‘ 二 ’ 一 ’ 一 或 ‘ · 若 以活的变化率为待求变量并将式 代人式 则可得 〔 〕 ’ 〔 〕 ” 〔 〕 ’ 。 〔 〕 ” 十 〔 〕 ’ 一 〔 〕 ’ 谧 一 〔 〕 ’ 〔 〕 ‘ 一 ’ 对于烤和温度而言 ， 它们的变化率存在如下 关系 若无相变发生 ， 则 〔 〕 〔 〕 〔 〕 因而有 〔 二〕 ’ 十 ‘ 〔 〕 ” ’ 〔 〕 ’ 其 中 ， 〔 〕和〔 〕 ” 是比热矩阵 。 若有相 变发生 ， 对于单元物质的相 变向题来说 即无 二相 区 ， 由于正在发生 相变 的节 点其温度不 发生变化 ， 故 、 两式不 再适用 。 但若采 用下述方法处理相 变节点所在单 元的热 导阵〔兀〕 ， 仍 可使两式成立 ， 例如 ， 如果 某单元的第 个节点处于相变 ， 则令该单元 〔尤〕 中第 行和第 列 中的元素为零 ， “ 二 。 意味 着可 以维持 节点的温 度不 变 ， 因 此 ， 将 关系式 代人方程 可得 〔 〕 △ 〔 〕 ， 〔 〕 〔 二 〕 ’ 〔 〕 ’ 一 〔 〕 ’ 一 〔 〕 ” 〔 〕 ’ 一 ’ 利用方程式 ， 本文在实 际计算 中采 用 ” 一 ” 一 的方法 求解方程式 可得〔 ，， 〕 ” ’ ， 然后用式 求 出〔 〕 ” 于 ’ ， 再 由治和温度的关系 曲线 可以 求 出〔 〕 ’ 十 工。 由上面的 推导可 以看出 ， 引 入三 时 间格式 以 后 ， 可以不需要迭代 而达到计算一次完 成的 效果 。 数 值 计 算 为了说 明上述数值计算方法的可行性 ， 本 文就一维半无限 空 间的凝 固问题和 二维无限大 方形角区 液体的凝固问题进行 了计算 ， 在计算 中 ， 采用双线性四边形等参单元 ， 单元热 导阵 的计算采 用补偿后 的 一点高斯积分 ‘ “ ’ ’ ， 同时 ， 还采 用 了集 中质量法处理 质 量 矩 阵 的 方 法 。 图 所示是本文采用 的一维半无限 空 间的简化模型 和有限元 网格划 分图 。 为 了 便 于 比 较 ， 本文采 用 如下 的物性数据 ， 即
X=0 X=L 图1一维离散形式 Fig.1 One-dimensional discretization form ks=2.22W/m·℃，k1=0.556W/m·℃ (pC)s=1.762x10J/m3.℃ (pC)1=4.226×108J/m8.℃ pL=3.38×108J/m3 同时，采用的计算区域和边界条件为， T,=10℃，Tw=-20℃，L=1m。 公式中采用的比例系数a=0.5。 图2所示就是计算结果与精确解的比较图，由图可以看出，二相交界面位置的计算结果 与精确解符合的很好。 图8所示是对二维无限大方形角区液体的凝固问题的数值计算结果与分析解之间的比 较，其计算条件为：初始时液体的温度为0且是液体的熔点，突然将边界温度降至一1，采 用的热物性数据为： ks=k1=1;Cs=C1=1 数值计算分别对L=1.561和工=0,25二种情况进行了计算，从计算结果看，也较好 地与分析解相吻合。 2.0 7:=283k Tw-293K Exacting solution o Caleulating result 1.5 0.2 1.0 。TL=0.25mT0=0.3℃ 0.1 Exueting soiution 0.5 oCaloulating resul L=1.5613mT0=0.0℃ 10 20 0 0.0 1x10's 0.0 0.51.0 1.52.0 xA 4at 图2本文钻果与精确解比较 图3本文结果与精确解比较 Fig.2 Comparison between the exacting Fig.3 Comparison between the exacting solution and the calculating result solution and calculating result 231
万 二 万 二 无 图 一维离散形 式 一 掩。 。 ℃ ， 一 ℃ 。 二 “ “ 一 ℃ 。 “ “ ℃ 户 同 时 ， 采用 的计算区域和边界条件为 ， ‘ ℃ ， 二 一 ℃ ， 。 公 式 中采 用 的比例系数 。 图 所示就是计算结果与精确解的 比较 图 ， 由图可 以看出 ， 二相 交界面位置的计算结果 与精确解符合的很好 。 图 所示是对二维 无限大方形角区 液体的凝 固问题的数值计算结果与 分析解 之 间 的 比 较 ， 其计算条件为 初始时液体的温度为 。 且是 液体的熔点 ， 突 然将边界温度降 至 一 ， 采 用的热物性数据为 。 一 。 一 数值计算分别对 二 和 二 。 二种情况进行了 计算 ， 从计算结果 看 ， 也 较 好 地与 分析解相吻合 。 陌旧忘 气 「一 一 八 刃三 ‘ 鬓 屯 一 奋 ’ 卜寸 、 ‘ 卜户洲 尸丫 夕尸 二 住 ，。 一 、 之互 万飞 ，飞资 一 丫 乓 长 、 冬 。 。 。 二 ， 丫 ， 砂 “ 号 ‘ 协 二 ” · “ 飞 二 副 … 图 本文 结果与精确解比较 几 立 零 ‘ 扩诬丽 图 本文结果 与情确解比较 几
4结 论 在将三时间格式的离散化方法引入有限元半离散方程的基础上，提出了一种求解烙法方 程的有效算法，从而取消了迭代，大大节省了计算量，使计算时间大大减少，从一维半无限 大凝固问题和二维无限大方形角区凝固的二个算例分析看，运算时间可缩短到1/15左右， (PC微机上实行，一维问题取16个节点，二维问题取144个节点)。而且可以看到，计算结 果能较好地与精确解的结果相吻合，可见本文提出的有效算法可靠而且实用。 参考文献 1 Shamsundar N,Sparrow.J.Heat Transfer,1975,8:333 2 Dalhuijsen A J,Segal A.Int.J.Num.Meth.Eng.,1986,23:1807 3 Liu W K,Belytschko T.Int.J.Num.Meth.Eng.,1984,20:931 4 Liu W K,Ong J S,Uras R Z.Comp.Meth.Appl.Meth.Eng.,1985, 53:13 邢钢炉外喷粉增硅脱硫技术 本系统采用了先进的微孔流化床技术，保证了喷粉工艺的稳定；另外，根据我国的国情 和该厂的实际情况，大胆采用未除湿的压缩空气作为喷粉载气，不仅大大地节省了生产成 本，同时在我国北方气候条件下，也保证了喷粉工艺要求，是我国在高炉铁水喷粉处理中使 用压缩空气载气最早和最好的。另外，当时炼铁和铸造行业小少人对喷粉增硅生铁的性能和 生铁质量持怀疑态度，针对这种担心，我们开展了对喷粉增硅生铁性能研究，结果发现在成 分相当的条件下，未经喷粉处理的高炉生铁的抗拉强度仅80~140MPa,而经过喷粉增硅生 铁的抗拉强度达到180~240MPa。这一研究结果，为高炉铁水的直接应用提供了依据。 两年来已喷粉增硅处理铁水7×104t,命中率86%，喷粉脱硫处理铁水7.2×104t,命中率 92,5%,生铁合格率由93,49%提高到99.85%，铸造生铁优级品率由85.18%提高到94.92%。 由喷粉增硅和脱疏所获净经济效益295万元。 本项研究成果在国内具有先进水平，而且设备结构简单适于中小厂家推广使用，投资省， 操作方便，有明显的经济效益、 232
结 论 在将三 时间格式的离散 化方法引入有限元半离散方程 的基础 上 ， 提出 了一种求 解给法方 程的 有效算法 ， 从而取 消了迭代 ， 大大节省了计算量 ， 使计算时间大大减少 ， 从一维半无限 大凝 固向题和二维无限 大方形角区凝 固的二个算例分析看 ， 运算时 间可缩短 到 左 右 ， 微机上实行 ， 一维 问题取 个节点 ， 二维 问题取 个节点 。 而且可 以看到 ， 计算结 果能较好地与精确 解的 结果相 吻合 ， 可见 本文 提 出的有效算 法 可靠 而且实用 。 参 考 文 献 “ ， ， ， ， ， ， ， ， ， ， ， 。 · ， ， 邢钢炉外喷粉增硅脱硫技术 本系 统采用 了先进 的微孔流化床技术 ， 保证 了喷粉工艺的 稳定 另外 ， 根据我 国的 国情 和该厂的实际情况 ， 大胆采 用未除湿的 压 缩 空 气作 为喷粉载气 ， 不仅 大大地节省 了生 产成 本 ， 同时在我 国北方气候条件下 ， 也保证 了喷粉工艺要求 ， 是我 国在高炉 铁水喷粉处理 中使 用 压缩空气载气最早和最好的 。 另外 ， 当时炼 铁和铸造行业不 少人对喷粉增硅生 铁的性能和 生 铁质量持怀疑态度 ， 针 对这种担心 ， 我们开展 了对喷粉增硅生 铁性能研究 ， 结果发现在成 分相 当的 条件下 ， 未经喷粉处理 的高炉生 铁的抗 拉强 度仅 ， 而 经过喷 粉 增硅 生 铁的抗 拉强度达到 一 。 这一研究结果 ， 为高炉 铁水的直接应用 提 供了依据 。 两年来 已喷粉增硅处理铁分仃 义 ， 命 中率 ， 喷粉脱硫处理 铁水 ‘ ，命 中率 ，生 铁合格率 由 纬提高到 。 ， 铸造生 铁优 级品率 由 。 提高到 肠 。 由喷粉增硅和脱硫所获净经济效益 万元 。 本 项研究成果在国内具有先进水平 ， 而且设各结构简单适于 中小 厂 家推广使 用 ， 投资省 ， 操作 方便 ， 有明显的经济 效益