多元多相体系化学平衡组成的计算

D0I:10.13374/j.issn1001-053x.1986.s2.007 1986年9月 北京钢铁学院学报 Special issue Journal of Beijing University 专辑2 of Iron and Steel Technology No2,1986.9 多元多相体系化学平衡组成的计算 何英英 张捷宇 程述武 摘要 本文从用计算机代替人工计算多元多相复杂体系化学平衡问题出发，对应用 计算机的Brinkley法，Vcs法及W hite法作」讨论和比较。指出White法在计算 上实现解决计算多元多相体系化学平衡问题的优点。并介绍了用White法编制的 计算多元多相体系化学反应平衡的计算机应用程序的结构及基本工作原理。列出等 题。其计算结果与位于联邦德国啊亨大学的治金热力学数据库、中国科学院化工冶 金研究所的无机热化学数据库的计算结果相一致。 The Calculation of Equilibrium Composition in Multicompon- ent and Multiphase Chemical Reaction System He Yingying Zhang Jieyu Cheng Shuwu Abstract In this paper we discussed Brinkley's method,VCS's method and White's method for the calculation of multicomponent and multiphase chemical equilibr. ium.White's method has the advantage of high calculation speed.It is no use fo White/s method determining the independent reactions and to sort the indepen- dent variables.The initial values of composition are selected arbitrarily.On th basis of these reasons mentioned above White's method is chosen in our compu ter program.In the case of treating multiphase system,White's method develo- ped by W.R.Smith is chosen to avoid the difficulty of numercial singularities when solving a set of linear cquations.The paper introduced construction and principle of the program.The calculated results are conformed with the results that calculated by the database of Theoretical Metallurgical Institute of Rhein- isch Westfalische Techisvhe Hochschule Archen West Germany and the Inorganic Thermochemistry Database of the Institute of Chemical Metallurgy Academia Sinica. ·42·

年 。 月 北 京 钢 铁 学 院 学 报 那 认 ‘ 专辑 ， ， 多元多相体系化学平衡组成的计算 何英英 张捷宇 摘 要 程述武 本文从用 计算机代替人工 计算多 元 多相复杂体 系化学 平衡问题 出发 ， 对应用 于 计算机 的 址 法 ， 法及 法作 讨论和比 较 。 指 出 法在 计算机 上 实现解 决计算 多 元 多相休 系化 学平 衡 问 题 的 优 点 。 并 介绍 了用 法编 制的 计算 多元 多 相休 系化学反应平衡的 计算机应用 程序的 结 构及基本工 作原 理 。 列 出 考 题 。 其计算结果与位于联邦德国阿亨大学的冶金热力学数据库 、 中国科学院化工 冶 金研究所的无 机热化学数据库的 计算结果相一致 。 ， 犷 ， 产 ‘ 王 ‘ · ‘ · ， 时 一 丁 竺 址 七 产 丫 ， ， 祖 。 。 。 了 ‘ 萝 时 。 DOI ：10．13374／j ．issn1001－053x．1986．s2．007

在科学研究和生产实际中经常涉及到计算不同条件下多元多相体系化学平衡组成的 问题。由于当前各学科的相互渗透，在解决这类问题中，代数方法得到了广泛的应用。 因为平衡混合物组成实质上是化学反应化学计量的函数，而这些函数应该能够使体系的 任何一个热力学势达到极值。求极值的每一个步骤都可以归结为线性代数的问题。这为 计算机解决多元多相体系化学平衡组成的计算，提供了理论根据。 对于一个简单体系来说，借助于用热力学求出平衡条件，可以手工计算其化学平衡 问题，但对于复杂体系来说，用热力学求出平衡条件计算其化学平衡问题所导出的超越 方程组的求解非常困难，手工计算往往不能胜任。 用计算机解决代数问题不仅有成功的经验，而且有成熟可靠的算法。从目前看计算 机解决化学平衡问题较通用的方法主要有：Brinkley16法，VcS4·7法及 White2'3'6]法三种。 1主要方法的比较 考查以下三种主要方法： 1.1 Brinkley法（又称平衡常数法）： 设体系中可能存在C个组分，将C个组分间可能发生的反应归纳为r个独立反应： ∑VB:=0j=1,2,…,r (1) i=1 平衡时 c c ZUiG1° i=1 RT—）=K, I ai-ii=expi=1 (2) j=1,2,…,r 以（2)式为基础，由反应平衡条件和物料平衡关系得到r+s阶联立方程： c i=1,2,…,C Ⅱa;uiH=ki (3) i=1 j=1,2,",r C ∑an=b1 1=1,2,…,S i=1 (1)、(2)、(3)式中，n为摩尔数，G为自由解，，B为分子式，a为活度，k 为平衡常数，U:表示j反应方程式中i组分的系数。对反应物v:为正，对生成物v为负。 a为i分子中1元素的原子数，b,为原料中l元素的总摩尔数。s为体系中的元素数。 联立方程组(B)可用Newton-一Raphson法求解。 1.2VCS法（又称反应进程法）： 将体系中所有可能发生的化学反应归纳为一组规范化的独立反应，再根据最小自由 能原理，同时求出每个反应的合理进程，使体系逐步向平衡靠拢。 恒温恒压下平衡时 ·43·

在科学研究和生产实际 中经常涉及到计算不 同条件下多元多相体系化学平衡组成的 问题 。 由于 当前各学科的相互渗透 ， 在解决 这 类间题 中 ， 代数方法得 到 了广泛的应用 。 扩 因为平衡混合 物组 成实质上是化学反应化学计量的 函数 ， 而这些 函数 应该能够使体系的 任何一个热力学势达到极值 。 求 极值的每一个步骤都可 以归结为 线性代数的问题 。 这为 计算机解决多元多相体系化学平衡组成的计算 ， 提供 了理论根据 。 对 于一 个简单体系来说 ， 借助于 用热力学求 出平衡条件 ， 可 以手工计算其化学平衡 问题 ， 但对于复杂体系来说 ， 用热力学求 出平衡 条件计算其化学平衡 问题所 导 出的超越 方程组的求解非常 困难 ， 手工计算往 往不能胜任 。 用计算机解决代 数 问题不仅 有成功的经 验 ， 而且有成熟可靠的算法 。 从 目前看计算 机解决 化学平衡 问题 较通 用的方法主要 有 厂 “ “ 、 法 ， 二 少法及 二 ’ ’ “ 〕 法三 种 。 主要方法的比较 考查 以下 三种主要 方 法 。 法 又称平衡常数 法 设体系 中可能存在 个组分 ， 将 个组分 间可能发生的反 应归纳为 个独 立反 应 名 ， ， ， 二 ， 平衡时 ” 二 ’ 一 ” ’‘ “ 二 一 ， ‘ ， 以 式为基础 ， 、 匆、 汀 二 ， ， … ， 由反应平衡条件和 物料 平衡关 系得到 十 阶联立方程 王二 ， ， … ， 。 一 “ ， ， … ， 乙 二 ， ， … ， 、 为平衡常数 ， 式 中 ， 为摩 尔数 ， 为 自由解 ， · 为分子式 ， ‘ 为活 度 ， 叭。表示 反 应方程式 中 组分 的系数 。 对反 应物 ，为正 ， 对生成物场 ，为 负 。 为 分子 中 元素的 原子数 ， ，为原料 中 元素的 总摩尔数 。 为体系 中的 元素数 。 联立方程组 可用 一 法求解 。 法 又称 反应进程法 将体系 中所 有可能发生的 化学反应归纳为一组规范化的独立反应 ， 再根据 最小 自由 能 原理 ， 同时求 出每个反 应的合理进程 ， 使体系逐 步向平衡靠拢 。 恒温 恒压下 平衡时

MinG=ni (Gi ( (4) i=1 为一组规范化独立反应的反应进程。5与关系为： n1=n:0-∑vi51≥0i=1,2,…,s (5) j=1 ni'=n°+5，≥0 j=1,2,…,r (6) 式(4)、(5)与（6)构成了一个完整的非线性规划问题，可用近似法求解。 1,3 White[°、1法（又称最小自由能法） 设体系中可能存在c个组分，其中第1,2，…m个为气相和溶液相物质，第m+1, m+2,℃个为凝聚相物质。在此我们假设每个凝聚相物质都是纯物质，则体系的总 自由能为： m c G= ∑ni(G:°+RTIna)+Σn,G;”(7) i=1 i-m+1 式中，a1=i为物质在所处相中的活度。 n r:为i物质在所在相中的话度系数。 n为i物质所在相中所有组分的摩尔数总和。 对于每一种元素都可写出一个物料平衡方程 ain=b11=1,2…,5 i=1 b是元素1在原料中的总摩尔数。在一个封闭体系中，无论发生什么化学变化或相 变化，体系中任一元素的总摩尔数不变。 在恒温恒压条件下，体系处于平衡的条件是：在满足物料平衡方程的前提下，使体 系的自由能G值为最小。这样，就把一个求多元多相体系平衡组成的问题，归结成一个 求有约束条件的最小化问题。 需要指出的是，在（7)式函数G的表达式中，不是所有变量都是独立的，因为这 些变量都已经由物料平衡方程联系起来，为找出函数G的极值，本应该首先借助体系中 相应的独立反应式，消去非独立变量后，再对独立变量求导数并令其等于零，找出极 值。但是，对于White法来说，它的便利之处恰恰就在于它不需要分出独立反应、独立 变量，不需要进行消去步骤，写出体系中的反应方程式。实质上，这些反应已经包括在 物料平衡方程中了。 使用Lagrange乘因子x'将有约束条件的最小化问题转化为无约束条件的最小化问 题。 Q=G-z1'（∑an:-b1) (9) 1=1 i=1 由(9)式可看出，π，'不仅仅是个数，而且要具有能量单位13)。其性质可由下 式得出： ·44·

， 君 烤 若为一组规范化独立反 应的反应进程 。 咨与 关 系为 二 “ 一 习 乙，》 ， ， … ， 产 “ 乙 ， ， … ， 式 、 与 构成 了一个完整的非 线性 规划问题 ， 可用近 似法求解 。 〔 、 ‘ 。 〕 法 又称最小 自由能法 设体 系中可能存 在 个组分 ， 其 中第 ， ， 一 个 为气相和溶液相物质 ， 第 十 ， 十 ， 一 个为凝聚 相物质 。 在 此我 们假设每 个凝聚 相物质 都是纯 物质 ， 则 体系的总 自由能为 名 ” 乙 ， ” · 式 中 “ ‘ 二 ，。 · 令 为 物质‘在所 处相 中的活度 。 ，为 物质在所 在 相 中的活度 系 数 。 为 物质所 在相 中所 有组分 的摩 尔数总和 。 对于每一种元素都可写 出一个物料 平衡方程 乙 一 ， … … ， 是元素 在原料 中的总摩 尔数 。 在一个封闭体系 中 ， 无论发生什么 化 学 变化或 相 变化 ， 体系 中任一元素 的总摩 尔数不 变 。 在恒温恒压 条件下 ， 体系处于平衡 的条件是 在满足 物料 平衡方 程 的 前提 下 ， 使体 系的 自由能 值为最小 。 这样 ， 就 把一个 求 多元多 相体系平衡组成的 问题 ， 归结成一个 求有约束条件的 最小 化问题 。 需要 指 出的是 ， 在 式函数 的表 达式 中 ， 不是所有变量都是独 立 的 ， 因为这 些 变量都已经 由物料平衡方程联系起 来 ， 为找 出函 数 的极值 ， 本应该首先借助体系 中 相应的独立反应式 ， 消去非独立 变量后 ， 再 对独 立变量 求导 数并令 其 等 于 零 ， 找 出极 值 。 但是 ， 对于 法 来说 ， 它的 便利 之 处恰恰就在 于它不需 要 分 出独 立反应 、 独立 变量石 不需要进行消去步骤 ， 写 出体系 中的反 应方程 式 。 实 质上 ， 这些反应 已经包括在 物料平衡方程 中了 。 使用 乘 因子二， 将有约 束条件的 最小 化问题转化为无 约束条 件的 最小化问 题 。 一 万 艺 一 ， 由 式可看出 ， 二 不仅仅是个数 ， 而且要 具 有能量 单位 〕 。 其性质 可 由下 式得 出 · ·

x,'=∑a1(G:°+RT1na)1=1,2,s(9a) i=1 式中，(a)c"s为由体系中C个化合物线性表示S个元素的系数矩阵。 由上式可知，1'为1摩尔1元素对体系的自由能贡献。（式(9a)的推导见附录）。 若物料平衡方程得到满足，则G=Q,Q最小时G也最小。当Q最小时需满足： 30÷0 i=1,2,C oni (10) 0=0 1=1,2,…S 2π1 将G在n:°(i=1,2,…)点展开为二级Taylor级数。 G 1 cc G-G( (ni-n:°)+ ni=nio i=1j=1 (n:-n,)(;-n1) 2G (11) on:ni n=n;° G的导数值如下， 1≤i≤m G 3n; =G:°+RTInai 32G onji =RT (Iny)+RT (1-1) oni ni n 22G -R olnyi) onionj j=n+1,m+2,…c(12) cni G (n j=1,2,“,m 2n3n, lnY）-RT是 0 m+1≤i≤c G =G1° 92G 9nI enn =0j=1,2,c 将方程(12)代入(11)，即可得G的二级Taylor级数展开式，将此展开式代人 (9)式并求导，则可求得： G0-RT∑xa1=0 i=m+1,…,c(13) 0= 1=1 oni +ina+ RT n0 i12.m() 式中：=1'/RT ·45·

万 ‘ 艺 ， 。 。 。 二 ， ， … 式中 ， · 为 由体系 中 个化合 物线性表示 个元素的 系数华阵 。 由上式可知 ， 二 产 为 靡尔 元素对体系的 自由能贡献 。 式 的推导见附录 。 若物料平衡方程得到满足 ， 则 二 ， 最小时 也 最小 。 当 最小时需满 足 叹 ， ， … 。 。 二 。 ， ， … 气 。 万 将 在 ，。 二 ， ， … 乃 点展开为二 级 尹 级数 。 止” 一 二 几 。 二 ， 一 。 “ 一 、 “ 、 二 。 。 ， 。 一 ‘ ” 名 · 名 的导数值 如 下 抓 铆镇 二 ，“ 。 一 、 拓 口 。 。 筱名 二 ‘舞 ， 口 ‘ ， 。 ， 丝型、 一 奋 ， ， … 。 一 一 胜 一 “ ， ， 一 ， 宁 二 。 · 。 立 。 么 。 ， “ “ ， ， ” 、 斗、 将方程 代人 ， 即可得 的二级 级数展开式 ， 将此展开式代 入 。 式并求导 ， 则可 求得 忿 “ 一 习 为 ， 一 ， 声 ‘了 一 一 公 ” 二 ， ， …皿 一 二 一砂 一 一生 、 式 中‘ 、 畜 了 ·

那“；8amb1=01=12…5（a5》由(14)式解出 =-ac+na+a。+n(名，a） 1=1 :、i=1,2,…,m :(i6) 对i=1,2,…,m求和有： m s m En.(=n RT +lna;)(17a) i=11=1 i=1 再将(16)式对溶液相各物质求和： m m n(Sa)三，n时(RT+ (17b) i=m'+11=1 i=m'+1 将(13)式两边乘上n:”并对i=m+1,…,c求和，然后与(17)式合并得 C s ni +Ina)+a (18) i=11=1i=1 RT i=5+1 xb,=已n°(gT+lna,）+.,ni°”(1 1=1 i=1 i=5+1 RT 将(16)式代入(15)式，则有 m an,9元a+,2,,a1n+n2,a1(D0)C i出1 j=1i=m+1 i=1 m =b+.a:n(+lna,）1=1,2,…, 1=1 m 设Y1=∑a1·ai·n:,经整理有： i=1 m m .Y1x1+zan+（-n0-1).∑an° j=1 i=m+1 i=1 .m =b-,an:+.∑ain°( (21)1=1,2,…,5 i=1i=1 RT +Ina 联立方程(13)、(17)、(19)、与(21)得到s+Φ级的线性方程组 (中为体系中的总相数)。 对方程组求解后，根据关系式 m=nc骨+名n(器+lna) 1±1； 求得n.再由下式 n1'=n°+A(n:-n:°) 求得下次选代的初值，进行下一次选代，直至满足精度要求。 ·46·

一 一 亡 娜 一 渝 “ 全 ， ，’ ” ’ 一 二 二 ， ， … ， 一 。 。 。 〔里华 、 ， 〕 。， 二， 一 十 。 。 找 一 ， 由 式解 出 艺 万 泊 二 ， ， ， 一 ， 玉 对 ， ， 一 ， 求 和 有 ‘ 万 万 军 一 、 ， ， 二 名 ，” 再 将 式 对溶 液 相 各物质 求 和 宜 丫 乙 才 ，决，， 二 ， · 盆谷 ‘ 一 ， 将 式两边 乘 上 、 ” 岁手对 二 ， … ， 求和 ， 然后 与 式合 个得 ” 艺 万 一二 名 二 ， 。 ，一 气 ， 豆了 名 “ 艺 二 兀 一 “ 。 一卜 、 ” ︸︸ 。 将 式 代人 式 ， 则 有 公 。 “ 二 兄 、兀 “ 一 ， 一 乙 影 一 ‘ ， ·沐 ︸ ” 一 ” 。 ， ， … ， 设 、 名 一 、 ， · 。 ， 经 整 理 有 艺 下，，万 ， “ 艺 ，， ， 二 十 一 一 ” 艺 十 】 二 一 ‘ 艺 联立 方程 、 一 一 一 ’ ‘影 · ，一， ， ， 一 ， 、 、 与 得 到 十 。 级的 线性方 程 组 。 为 体系 中的总 相数 。 一 “ ” 二几 一 求 得 ， 再 由 一下式 对方程组求解后 ， 根据 关系式 万 节 ， 〕 ，尸 二 ，“ 入 一 “ 隶 得下次 迭 代的 初值 ， 进 行下一次 迭 代 ， 直至 满 足精度要 求 。 味

在这里有一个问题要注意，若我们选择的初值n°(i=1,2…,©)精确地满足物料平 衡方程，就可以在得到方程组的梢确解后，知道下次迭代值也能精确地满足物料平衡方 程。实际上，一般不易于获得方程组的精确解，计算机的机器误差限制了解的精度。因 此，即使：°满足物料平衡方程，1‘也不能服从物料平衡方程。经过多次的迭代，机器 误差的累积，即将导致White法无法解决平衡计算问题。为避免产生这样的情况，我们 对方程组作了如下处理。 由物料平衡方程可导出 ∑at(n;-n,)-c1=01=1,2,…,s i=1 式中，c=a,n,°-b1 i=1 将(16)式代入上式中，就可以得到一个方程组： m c m 2yx,+2an:+(--1)∑an j=1 i=m+1 n i÷1 m m =bi-,,ain,°+∑ann,( RT +1na:)-c(22) i=1 i=1 1=1,2,,5 由上式看出，c1在精确地符合物料方程时，恒等于零，方程（22)就化为方程(21)。 而当°、n;'不能很好地符合物料平衡方程时，c项就包括了由此产生的误差。在方程 的解趋于正确解时，c,→0，由于©的存在，就不必要求每次送代的初值都满足物料平衡 方程。 White法与Brinkley法和Vcs法相比，具有不必找出独立反应，划分独立变量的优 点，并且初值的选取对计算影响不大。鉴于White法以上优点，我们选用White法作为 求解多元多相化学反应平衡组成的基本算法。用FORTRAN语言编制了多元多相反应体 系平衡组成的通用计算程序。 2多元多相体系化学平衡组成计算程序 在热力学数据库应用系统中，存在着若干个应用程序，人们可以根据自己的工作意 图投入不同的命令，使用不同的应用程序，并通过应用程序的执行得到计算结果。本程 序就是应用程序之一。它的功能是找出平衡化合物，计算多元多相化学反应的平衡组 成。 为了便于对程序的维护、扩充和更改，此程序与其它几个应用程序是独立的，没有 互相调用的关系。 本程序由21个子程序组成，大约四千条程序。它受控于监控程序，并要通过系统程 序从数据文件中读取热力学数据。程序的主体结构如下： ·47·

在这里 有一个 问题要注意 ， 若我们选择的初值 拓 。 二 ， … 精 确地满 足物料 平 衡方程 ， 就可以在得到方程组的精确解后 ， 知道下次迭代值也能精 确地满足物料 平衡方 程 。 实际上 ， 一般不 易于获得方程组的精确解 ， 计算机的 机器 误差限制 了解的精度 。 因 此 ， 即使 。 ” 满足物料平衡方程 ， 产 也不能服从物料平衡方程 。 经 过 多次的 迭 代 ， 机器 误差的累积 ， 即将导致 法无 法解决平衡计算问题 。 为避 免产生这 样 的情况 ， 我们 对方程组作 了如下处理 。 由物料 平衡方程可导 出 乞 一 一 一 二 ， ， … ， 式 中 一 艺 边 ” 式代 入 上式 中 ， 一 将 就可 以得 到一个 方程组 兀 。 一 一 ” 三一 。 乙 一 ” 二 艺 ，一 。 二 箭 一 ‘ 一 ，一 ‘ ， 由毕 看 出 ， “ 在 ” 云。 精确 地符合 物料方程时 ， 砂 而 当 尹 、 厂不能很好地符合物料 平衡 方程时 ， ， ， … ， 恒等于零 ， 方程 就 化为方程 。 ，项就包括 了由此 产生 的 误 差 。 在 方程 的解趋于正确解时 ， ，， ， 由于 的存在 ， 就 不必要 求每 次迭 代的初值都满 足 物料平衡 方程 。 法与 法和 法 相比 ， 具 有 不 必 找 出独 立反应 ， 划 分独立变量的 优 点 ， 并且 初值的选 取对计算影 响不 大 。 鉴于 法以上优 点 ， 我们 选 用 。 法作为 求解 多零多相化学反应平衡组成的基本算法 。 用 语言编制 了多元多相反应体 系平衡组 成的通 用计算程 序 。 一 多元多相体系化学平衡组成计算程序 ， 在热力学数据库应 用系统 中 ， 存 在着若干个应用程序 ， 人们可以根据 自己的工作 意 图 投入不 同的命 令 ， 使用不 同的 应用程 序 ， 并通 过应用程序的执行得到计算结果 。 本程 序就是应用程序之一 。 它的功能是找 出平衡 化合 物 ， 计算 多元 多 相 化学反应的 平衡组 成 。 为 了便于对程 序的维护 、 扩充和更改 ， 此程 序与其它几个应用程 序是纯立的 ， 没有 互 相调用的关系 。 本程序由 个子程 序组成 ， 大约四千条程序 。 它受控于监控程序 ， 并要通过系统程 序从数据文件 中读取热力学数据 。 程 序的主体结构如下

READ 查找包含系统元素的化合物子程序 EQUI HSG·· SGMIX SORT GASOL 计算子程序 OPUT 程序的工作原理如图1： Start EQUI 图1总流图 READ Fig.1 General flowchart Calculate? Y SGMIX End EQU1:为了用户使用上的方便，当用户需要找出数据库中某一系统元素组成的平衡化 合物时，可投入EQUIFIL命令。此命令不做计算，仅仅查找化合物。如果用户既要找 出化合物，又要计算多元多相化学反应的平衡组成时，投人EQUI OMP命令。EQUI程 序就是根据用户投入的命令，按照用户的意图决定程序的走向。 READ,查找化合物程序。 into from EQUI SGMIX:.多元多相化学反应的平衡组成的计算子程序。 Decade the syatem ·相n@teB 3· 各主要程序段的功能及工作原理 Find syatem compounds in file 1 3.1READ程序的功能是调用系统子程序在数据文件一， 数据文件二中检索该系统元素的平衡化合物。工作原理如图 Solution? 2。 Y Find ayatem compounde 1n'f11a2 Print syetem conpounds 图2READ程序流程图 …Delet7 Fig.2.Flowchart of READ Program Y Cal1 RDNUMS Return ·48

一 “ 奴” 妙 包含系统元素的 化合物子程序 一 计算子程 序 了 了 ， 程序的 工 作原理如 图 图 总流 图 。 为了用 户 使 用上的方便 ， 当用户需要 找 出数据库 中某一 系统元素组 成的平 衡化 合物时 ， 可投入 命令 。 此命令不 做计算 ， 仅仅查找化合物 。 如 果用 户 既要找 出化合 物 ， 又要 计算多元多相化学反 应的 平衡组成时 ， 投入 命令 。 耳程 序就是根据用 户 投入 的命令 ， 按照 用户的意 图决 定程 序的 走 向 。 查找 化合物程序 。 一 多元多相化学反 应的 平衡组成的计算子程序 。 · 各主要程序段的功能及工作原理 人 ‘ 性 程序的功能是 调 用系统子 程序在数据文 件一 ， 数据文件立 中检索该系统元素的乎衡化合物 。 工 作原理如 图 勺叮一了 严 塑竺甲 ‘尹士，二一， ， · ‘ ， 乙 厂 名

Into from EQUI 在执行本程序时，用户要在命令参数system中，提 PTOP=1 供所要计算的系统的全部元素，最多可达20个。该程序 NI0T=1,E+6 TLIM-20 首先将命令的system参数中的元素一个个的分解成纯物 Input T 质，并根据分解出的元素，顺序的从文件一中找出含系 统元素的所有化合物。如果各元素在某种溶剂中溶解时 Solution？ Call HS ,还要在文件2中顺序找出系统元素及此溶剂组成的所有 Cell SORT 化合物。将这些化合物的分子式及最高挡温度存于数组 Read the program 之中并打印出来供用户参考。用户可对不予考虑的化合 execution directive Coll.GASOL 物进行删除。删除功能由RDNUMS程序实现。当上述 工作完成之后，程序将化合物的分子式及它们的原子矩 Call OPUT 阵写人文件DAT供计算时使用。 End 3.2 SGMIX程序的功能是计算多元多相化学反应的平. 图3 SGMIX程序流程图 衡组成。它是通过调用几个辅助子程序和计算子程序实 Fg.3 Flowchart of SGMIX Program 现的。工作原理如图3。 程序首先将反应物质的总重要、气压等置初值。为了避免在某些情况下过分延长计 Input from SGMIX 算时间，程序在用户没有给出计算时间时，给出每次平 衡计算的最高时间限制。然后由用户输入反应温度值， Sort the system compounds into a gas phase and a condensed phase 系统中各物质将自动按此温度区分各个相区，如气相或 凝聚相。当有稀溶液存在时作为附加的另一相，·接下来 Get the date for the substences 执行HSG子程序。 Call GVALS HSG程序是SGMIX的一个计算子程序。它的功能是 计算化学反应中溶液的超额自山能△G严。 Solution？ Y SORT程序是SGMIX的一个辅助子程序。工作原理 Call FNDEPS 如图5。它的功能是将系统元素组成的化合物按相分类 End 排序存于数组之中，然后根据化合物在数组中的位置顺 图4SORT程序流程图 序从文件一或文件二中取出计算时所需的参数，组成参 Fig.4 Flowchart of SORT Program 数表。 GVALS程序是SORT程序的一个计算子程序。它的功能是根据参数表中提供的数据 计算温度T下系统中每个化合物的△G,CP,H,G,S。 FNDEPS程序的功能是当有稀溶液时。计算两物质相互作用系数。 在SGMIX程序中为了便于用户对温度、压力、物质总重等各参数进行修改或为计 算选定初始反应物和初始值，提供了十种不同的处理方式，用户可以输入1到10数字中 的一个来决定选用的处理方式。这些方式如下： (1)进一步删除系统化合物 由于用户指定的温度各有不同，因此各物质及其存在的相区也有区别。计算机在执 行此处理时先打出化合物清单（各化合物的序号、分子式及温度范围），供用户参考。 再次为用户提供删除机会。 (2)给定反应物质的初值 49·

图 程序 流程图 在执行本程序时 ， 用户要 在命令参数 中 ， 提 供所要 计算的系统的全部元素 ， 最 多可达 个 。 该程序 首先将命 令的 参数 中的元素一个个的分解成纯物 质 ， 并根 据分解 出的元素 ， 顺序的从文件一 中找 出含系 统元素 的所 有化合物 。 如果 各元素在某种溶剂 中溶解时 ，还要 在文件 中顺序找 出系统 元素及此溶 剂组成 的所有 化合 物 。 将这些 化合 物的分 子式 及最 高挡 温度存 于数组 之 中并打 印出来供用户参考 。 用户可对不 予考虑的 化合 物进 行删 除 。 删 除功 能 由 程序 实 现 。 当 上述 工 作完成 之后 ， 程序将化合物 的分 子式及 它们的原 子矩 阵写人文 件 供计算时使用 。 程序的功 能是计算 多元 多相化学反应的 平 衡组 成 。 它是通 过调 用几个辅助子程序和计算子程序实 现的 。 工作原理如 图 。 程 序首先将反应物质的总重要 、 气压等置 初值 。 为 了避免在某些情况下 过分 延长计 七 日 呻 日 吕 图 程序流程图 】 算时 间 ， 程序在用 户没 有给 出计算时 间时 ， 给 出每 次平 衡计算的 最高时 间限制 。 然后 由用 户 输人反 应温度值 ， 系统 中各物质将 自动 按此温度 区分 各个 相区 ， 如气相或 凝聚 相 。 当有稀溶液存在时作为附加的 另一相犷 接下 来 执行 子程序 。 程 序是 的一个计算子程序 。 它的功能是 计算化学反 应 中溶 液 的超额 自山能 △ 。 程序是 的一个辅助子程 序 。 工作 原 理 如图 。 它的功 能是 将 系统元素组成 的化合 物按 相分 类 排序存于 数组之 中 ， 然后根据化合物在数组 中的位置顺 序从文 件一或文 件二 中取出计算时所需的参数 ， 组成参 数表 。 程序是 程序的一个计算子程 序 。 它的功 能是根 据参 · 数表 中提供 的数据 计算温度 下系统 中每个 化合 物的 △ ， ， ， ， 。 程序的功能是 当有稀溶液时 。 计算 两物质相互作 用系数 。 在 程序中为 了便于用 户对温 度 、 压 力 、 物质总 重等各参数进 行 修 改或为 计 算选定初始反 应物和初始值 ， 提供 了十种不 司的处理方式 ， 用 户可 以输入 到 数字 中 的一个来决定选用的处理方式 。 这些方式 如下 进一步删除系统 化合 物 由于用户指定的温度各有不 同 ， 因此 各物质 及 其存在的 相区也有 区别 。 计算机在执 行此处理时先打 出化合物清单 各化合 物的序 一 号 、 分 子式及温度范围 ， 供 用户参考 。 再次 为用户提供删除机会 。 给定反应物质 的 初值

， 用户可以选择多达20种反应物的初值输入计算机，这包括初始反应物的数门，反应 物的分子式，反应物初始量。当这些初值量输入计算机之后就进入了多元多相化学反应 平衡组成计算子程序GASOL,这将在下面做介绍。 (3)可改变反应物质的总重 在用户不作改变时，总重的缺省值是10°克。 (4)用上次计算时的初始反应和初始量作为本次计算的初值。 用户可以改变状态（温度、压力等），但不改变反应物及初始量。 (5)可以改变一种初始反应物及初始量。 (6).改变反应温度. 系统中的各化合物将自动地按该温度下各物所的隐定重新分飞。 (7)改变总压 在用户不作改变时，总压的缺省值是一个大气压。 (8)给定反应物和初始量，并由用户输入估计平衡值。 (9)恢复程序原来对压力，物质总重等缺省值的指定。 (10)命令终让· 以上提供十种处理方式使用户可以对化学平衡化合物在不同的温度，不同的初始反 应物，不同的初始值及不同的压力进行计算，用户不需要从头运行本程序，这显示了木 程序的灵活性。 GASOL程序是SGMIX程序的一个重要的计算子程序。它的功能是计算化学反应的 平衡组成。比较三种主要算法时已经对White法作了详细的介绍。本程序是根据White 法编写的程序。因此它的基本算法请参见前述Whte法。这里仅对程序中的主要步骤进 行说明。 （1)凝聚相的选择： 文献6指出，用White法解决多元多相化学平衡问题时，遇到下面几种情况，其线 性方程组的系数矩阵发生奇异现象而导致计算失败。 ①在计算过程中某凝聚相趋于消灭时，线性方程组系统矩阵发生奇异现象而导致计 算失败。为解决此问题，采用选择凝聚相的方法，即若一相中的物质总量小于10~19 摩尔，则在系统矩阵的计算中去掉此相，而路掉的这一相物质的组成在方程组求解后， 由关系式 :=exp(-+2,an) (23) 求出其活度值。式(23)的推导见附录。若还未得到问题的正确解，则需求出略去相的 组成，进行下一次迭代。由于略去相的物质总量小于10-1·摩尔，可以认为活度系数 Y:=1,故有 n=nexp(8T+∑an +ai) (23a) 由式(23a)可求出略去相的物质组成。 ·50·

用 户 可 以选择 多达 种反应物的 初值 输人计算机 。 这包括 初始反应物的数月 ， 反应 物的分子式 ， 反应 物初始量 。 当这些 初值量输人计算机之后就进 入 了多元 多相化学反 应 平衡组成计算子程序 ， 这将在下面做 介绍 。 可改 变反 应物质的总重 ’ ‘ 一 ‘ 一 在 用户 不作改变 时 ， 总重的缺 省值是 “ 克 。 用 上次计算时的 初始 反 应和初始量作 为本次计算的 初值 舀 一 用 户可以改变状态 温 度 、 压 力等 ， · 但 不 改 变反应 物及 初始量 。 一 可 以 改 变一种初始反应 物及 初始量 。 改 变反 应温度 ‘ 乙 一 系统 中的 各化合 物将 自动地 按该温度 下 各物质 的稳 定 相 重新分 哭 。 改 变总压 一下 一 在用户 不作 改变时 ， 总 压 的缺省值是一个 大 气压 。 、 一 给定反 应 物和初始量 ， 并 由用户 输入估 计平衡俏 。 恢 复程 序原 来对压 力 ， 物质总 重 等缺 省 值的指定 。 命令终 止 · 以上提供十种处理 方式使用 户可 以 对 化 学 平衡 化合 物 在不 同的温 度 ， 不 同 的 初始反 应物 ， 不 同的 初始值及不 同的 压 力进行 计算 ， 用 户不需 要从头 运 行 本程 序 ， 这显示 子本 程序的灵 活性 。 一 ’ 卜 程序是 程 序的一个 重要 的 计 算子程 序， 它的功 能是计算化学反 应 的 平衡组成 。 比较三种 主要算 法时已经对 瓦 法作 了详细 的介绍 。 本程序是根 据 法编写的 程 序 。 因此 它的基本算 法请参见前述 法 。 这 里 仅 对程序 中的 主 要 步骤进 行 说 明 。 ‘ 凝聚相的选 择 文献 一 ，指 出 ， 用 切法解决 多元 多 相化学平衡 问题 时 ， 遇 到下面几种情 况 ， 其线 性方程组的 系数矩阵发生奇异现象而导致计算失败 。 ①在计算 过程 中某凝聚相趋于消灭时 ， 线性方程 组 系统矩阵发生奇异现象而导致 计 算 失败 。 为解决 此问题 ， 采用选 择凝聚 相的方 法 ， 即 若一 相 中 的物质 总量小 于 一 ‘ 。 摩 尔 ， 则在 系统矩阵的计 算中去掉此相 ， 而略掉 的这一相物质的组成在方程组求 解后 ， 由关 系式 一 旦二 兀 求 出其活度值 。 式 的推导 见附录 。 若还未得 到 问题的正 确解 ， 则需录出略去相的 组成 ， 进 行 下一 次迭 代 。 由于略 去 相的 物质 总 量 小 干 一 ， 。 摩 尔 ， 可 以认 为活度系 数 丫 ， 故有 ， 一 ‘ 盖 一下刃丙丙一 工吸 万 一 二 由式 可 求 出略去 相的 物质组 成

②线性方程组在有多于S-1个凝聚相物质参加运算时，其系数矩阵就发生奇异现 象。用下面的方法解决此问题，在开始计算线性：方程组的系数矩阵时，我们假设所有凝 聚相物质量为零、在方程组得到解以后，用下式 △G;=G:°/RT-∑a1i=s+1,,c 1=1 计算每个凝聚相物质的自由能。在下一次迭代方程组的计算过程中，考虑具有最小△G 的凝聚相，并重新计算平衡量。在新的平衡下对剩余的凝聚相物质的△G:值进行比较， 将最小者加入到新的迭代过程中，这个步骤进行下去直到所有剩余凝聚相物质的△G:≥ 0,计算出平衡组成。对于方程组中未考虑的凝聚相物质，其平衡量由式(23)计算。 (2)寻找合理的迭代步长因子入 选择迭代步长因子需满足下面两个条件： ①n;'=n:°+入(n:-n:)≥0 ②λ的选择使得：'不超过体系自由能的极小值。即体系的自由能在选代方向导数不 大于零。 d(CRDr+ing(t() dλ “z(n,+λ(n;'-n) 式中工表示组分所在相所有组分相应项之加。 我们选择满足以上二条件的最大的入作为步长因子，求得下次迭代初值：'： n;'=n:°+λ(n:-n.) 在求出平衡组成达到精度后，入=1。 ③检查计算是否收敛 设线性方程组的解为Y={y:}(i≠1,2，…s+中)，e,、eg为任意给定的精度要 求。 若 y二y.se 且 ;ni-n°≤ea ni 则说明计算出的结果为正确解，由计算机输出平衡组成。 (4)初值的选取 在方程组中由于©：(=1,2，…，s)的存在，使得初值的选取对计算的收敛性影响不 大，故可任意选取，只要n:°≥C。 (5)热力学模型1〕 在常压、，高温的条件下，气体按理想体系考虑，蟹聚相的个数按其凝聚化合物的个 数计算，熔体的超额自由能用Redlich-Kister-Muggianu模型计算，选纯物质为标准态。 Redlich-Kister-Muggianu模型， …51·

②线性 方程组在有 多于 一 个凝聚 相物质参加运算时 ， 其 系 数矩 阵就 发 生 奇 异现 象 。 用 下面 的方 法解决此问题 ， 在开 始计算 线性方程组 的系数矩阵时 ， 我们假设所有凝 聚相物质量 为零 ， 在方程 组得 到解以后 ， 用下式 △ 一 ‘万 ， … ， 计算每个 凝聚 相物质的 白山 能 。 在下一 次迭代方程组的计算 过程 中 ， 考虑具 有最小 △ 的凝 聚相 ， 并重新计算平衡 量 。 在新 的平衡下 对剩余的凝 聚相 物质的 △ 值进行比较 ， 将最小 者加入 到新的迭 代过程 中 ， 这个步骤进 行下去 直到所 有剩余凝聚 相物质的 八 ‘妻 。 ， 计算 出平衡 组 成 。 对 于方程组 中未考虑的凝 聚 相物质 ， 其平 衡 量 由式 计算 。 寻找合理 的迭 代步 长 因子入 选 择迭代步长 因子需满 足下面两个条件 ① 盛‘ ‘ 。 、 ’ 盛 一 ， 妻 。 ②入的选择使得 。 ’ 不超 过体系 自由能的极小值 。 即体系的 自由能在迭 代 方 向 导 数 不 大于零 。 入 一” 歹 ， 〔补 十 。 少 一上入典乓牛西互 一 〕 ， 一 五 ， 乙 人 吸 ‘ 一 〔 式 中艺 表示 组分 所在 相所有组分 相应项 之加和 。 我 们选 择满 足以上二 条件 的 最大 的入作为步 长因子 ， 求得下次迭代初值 ‘ ‘ ‘。 入 一 五 ‘ 在求 出平衡组 成达 到精度后 ， 入 。 ③检查计算是 否收敛 设线性方程组 的解 为 厦 “ ， ， “ · 十 中 ， 。 ， 、 。 为任 意给定 的精度 要 求 。 若 ‘ 一 ” 簇 一 玉 ￡ 则说明计算 出的结果 为正确解 ， 由计算机输 出平 衡组 成 。 初值的选 取 在方程组 中由于 二 ， ， … ， 的存 在 ， 使 得 初值的选 取 对计算 的收 敛性影响不 大 ， 故可任意选 取 ， 只要 “ 妻 。 热力学模型 ’ 〕 在常压 、 ‘ 高温的条件下 ， 气体按理想体系考虑 ， 凝聚相 的个 数按其凝聚 化合 物的个 数计算 ， 熔体 的超额 自由能用 一 州 。 卜 模型计算 ， 选纯 物质为标准态 。 ， 一 一 ‘ 模型