水和水蒸气密度的跨区域计算

中国热物理学会 工程热力学与能源利用 学术会议论文 编号：101088 水和水蒸气密度的跨区域快速计算方法 安宾1, 王智学1, 王晓东1*,段远源2 1 北京科技大学机械工程学院热能工程系，北京 100083 2 清华大学热科学与动力工程教育部重点实验室，北京 100084 010-62321277，wangxd99@gmail.com 摘要: 采用双二次样条插值的方法，建立了水和水蒸气密度的双二次样条方程ρS(p,T)。有效计算区域 为：温度 273.16 K到 1073.15 K；压力 0.001 MPa到 100 MPa，采用双网格技术样条模型，能达到很高 计算精度，计算速度比IAPWS-IF97 快 10-540 倍。本文对热力学面进行整体拟合，避免了分区导致的 区域边界数值不一致问题。 关键词：水和水蒸气；密度；样条插值；快速计算； 0 引言 水和水蒸气在工业生产特别是电厂热能动力工程以及相关研究领域是一种最普遍 最常用的介质，其热物性的精确、快速征直接关系到生产效率、经济性和安全性。尤 其电厂蒸汽动力循环的和优化以及运行监控，需要快速而精确的热物性计算方法。 随着科学的发展，对水和水蒸气性质的精确描述也在不断更新发展。由于国际水和水蒸 气协会(IAPWS)公布的“通用和科学用水和水蒸气性质计算公式”（IAPWS-95）计算极为 复杂，为满足工业应用需求，国际水和水蒸气协会(IAPWS)于1997年发布了 “IAPWS水 和水蒸气热力性质1997工业用公式”，简称为IAPWS-IF97[1]。IAPWS-IF97数据来源基于 IAPWS-95公式的计算结果，方程形式进行了适当简化。但是IAPWS-IF97的计算分为多 个区域，在边界处存在不一致性，同时在迭代计算时，由于计算公式被多次调用时，计 算速度慢、稳定性差，有时甚至得不到收敛解。 先进电力循环和过程的设计和优化需要快速和精确的热物性计算方法，因为工质热 物性需要在过程计算的内部迭代循环中频繁使用，有时甚至被反复调用多达 106或 107 次。表格查询法是实现快速计算的有效方法。例如典型的表格泰勒级数展开法（TTSE）， 计算速度比IAPWS-95 快 100 到 200 倍，且可达到较高计算精度[2][3]。尽管TTSE计算精 度和速度较高，正反函数具有较好的数值一致性，但其缺点是不能连续地表征整个热力 学面。 发展新的水和水蒸气热物性快速计算方法具有迫切的需求，新方法需要满足如下几 个条件：1、具有足够的计算精度；2、计算速度比 IAPWS-IF97 快；3、尽可能减少热 力学面的划分区域数目，有效解决各区域边界的一致性问题；4、可连续表征热物性的 整个热力学面。本文发展了基于样条函数的水和水蒸气热物性计算方法，发展了样条插 值模型。 1 样条插值模型 针对目前水和水蒸气快速计算方法存在的问题，二维样条模型是有效的解决途径之 一，其不仅具有表格查询法的优点，数据的预处理可节省大量计算时间，而且能够连续 的表征整个热力学面，同时具有高精度和高计算速度的特点。 基金项目：国家重点基础研究发展计划资助（No. 2009CB219803） 与一般的分段逼近方法相比，样条模型具有如下优点：1）样条函数及其导数在各 节点处连续；2）样条函数为低阶多项式，计算工作量小，速度快，而且有较高的灵活 性，在逼近复杂函数时也较少出现振荡现象；3）样条插值模型的系数可通过计算机预 处理获得并存储，能节省大量计算时间； 4）样条函数及其反函数计算具有完全的数值 一致性。 根据不同计算需求，可定义待求的物性参数z为另外两个参数x和y的函数。在x-y热 力学面上建立如图 1所示的网格，其中，m和n分别表示沿x和y方向的格点数目，整个热 力学面包括m×n个格点和(m-1)×(n-1)个子矩形[4]。 图 1 网格划分 本文采用矩形网格的优点在于算法简单，仅通过给定点坐标(x,y)即可确定其所在的 子矩形位置。本文选用的双二次多项式样条函数形式为： ∑∑ = = −− −−= 3 1 3 1 11 )()(),( k l l j k iijklij yyxxayxG (1) 式中，aijkl为待定的样条函数系数。显然每个子矩形包含 9 个待定系数，整个区域共有 9×(m-1)(n-1)个系数需要确定，因此必须给出同样数量的约束条件。 不同的约束条件可以构成不同的样条函数形式。为提高计算精度，引入了双网格计 算模型[5]，如图 2。 图 2 双网格模型网格划分 在格点 (i,j)处必须满足插值的要求。要保证 (∂G/∂x)y在 (i′-1,j)和 (i′,j)的导数以及 (∂G/∂y)x在(i,j′-1)和(i,j′)的导数连续。其次，要保证在四个角点(i′-1,j′-1), (i′-1,j′), (i′,j′-1) 和 (i′,j′)处，由样条函数Gij (x,y)和由原始方程计算得到的二阶交叉导数(∂G2/(∂x∂y))相等。通 过对边界节点导数连续所得到的方程组进行求解可确定样条多项式的系数。采用上述方 法确定了 9×(m-1)(n-1)个样条函数系数后，可将系数值存储起来，建立相应的数据库， 在计算热物性时，只要确定了计算点所处的网格位置，就可从数据库中调出相应的系数 aijkl，利用样条插值函数计算。网格的密度会对样条方程与标准方程之间的拟合精度有 影响，网格划分越密，拟合精度越高。而网格密度与计算时间无关，计算样条函数系数 为预处理过程，只需一次即可，因此大大提高了热物性的计算速度。 2 双二次样条函数 ρ(p,T) 密度与温度、压力的关系是流体最基本的热力学性质，使用状态方程法计算水和水 蒸气密度计算时一般需要多次迭代，为验证本方法的可行性，本文首先选取 p、T 作为 自变量，建立密度 ρ(p,T)的快速计算公式。ρ 为 p 和 T 为变量的函数，对于 p-T 平面上 给定了一个矩形区域: ( ){ }dTcbpayx ≤≤,≤≤:,=ρ (2) 在这个矩形区域上建立如图 3所示的双网格模型。整个矩形平面被均匀划分为m×n 个相同的子矩形，子矩形的四个角点为格点，每个子矩形的中心为节点，整个矩形区域 包括（m+1）×（n+1）个格点和m×n个节点。这些节点和格点便构成了两套网格，即节 点网格和格点网格。 图 3 ρS(p,T)函数子矩形示意图 在p-T热力学面上的计算区域内，通过IAPWS-95 计算公式来给出约束条件中网格数 据点的值及边界点的导数值[6-[7]。本文将双网格双二次样条模型应用于在p，T热力学面 上水和水蒸气的密度的计算，并通过与IAPWS-IF97 模型的对比分析样条函数的计算精 度和速度。 2.1 网格划分方法以及边界条件的确定 本文的有效计算区域为温度 273.16 K-1073.16 K，压力 0.001 MPa-100 MPa。将这个 区域划分为 100×100 的均匀网格进行样条插值试算，结果表明在不同区域拟合效果有差 异，在饱和线附近，计算偏差显著高于其它区域。因此在不同区域应采取不同的网格划 分密度，本文采用的划分区域如图 4所示。 图 4 样条函数ρS(p,T)全区域网格划分 在各个子区域内采用均匀网格，计算偏差较高的区域网格进行加密。对于每个计算 区域（子区域），在网格划分之后还需给出边界条件，每个子矩形都有 9 个待定系数， 必须给出 9 个约束条件才能求出 9 个待定系数。双二次样条模型没有必要给出所有子矩 形的 9 个约束条件，只需给出整个子区域的边界条件即可，包括：子区域边界四个角点 上待求热物性的关于 p、T 的二阶混合偏导值，左右边界各个子矩形边界中点的待求热 物性对 T 的一阶偏导，上下边界各个子矩形边界中点的待求热物性对 p 的一阶偏导。同 时给出的定解条件还有各个子矩形中心点即格点的待求热物性的值。本文中边界条件是 用精度为四阶的有限差分的方法给出的。 2.2 饱和线附近区域的处理 在p、T热力学面上，跨饱和线水和水蒸气的密度是不连续的。对包含饱和线的矩形 区域进行计算得到的样条函数偏差较大，而且越靠近饱和线偏差越大。因此，本文在饱 和线附近区域采用了非规则区域的网格划分方法，如图 4中的 5、6、9、10 等计算区域， 它们的非规则边界是饱和线。采用边界规范化的方法将这些不规则区域转换成计算平面 上的矩形区域，然后可将双二次样条模型直接应用于该区域，基本原理如图 5。 图 5（a）是一个不规则的区域，右边界为二次曲线，采用下列坐标变换把右边界规范 化： )(,, 1maxmax yfxxxy −=== ξη (3) 这里的yt为右边界节点的y值，变换后的区域形状如图 5（b）。对于类似的非规则区 域，只要给出了不规则边界上y与x的关系式，就可以用这种方法进行变换，得到规则的 矩形区域。 （a） （b） 图 5 边界规范化示意图 图 4中的 5、6、9、10 等以饱和线为非规则边界的计算区域，可以用上面介绍的边 界规范化方法进行处理，非规则边界的方程即饱和区温度方程Ts（p）。 2.2.1 近临界区的处理 在压力 22.064 MPa 到 30 MPa 的近临界区，密度的变化剧烈。对该区域直接进行样 条插值拟合，计算偏差较大。尽管跨拟临界温度线水和水蒸气的热物性ρ的值都是连续 的，但在其附近三个热物性的一阶偏导明显增大即热物性的值变化剧烈，近临界区的这 一特性类似于饱和线附近区域。 近临界区内通过网格加密的方法仍然不能很好的解决样条函数偏差较大的问题，本 文采用了和饱和线附近区域同样的处理方法，以拟临界温度线为非规则边界，应用边界 规范化的方法将非规则区域变换为规则区域，再应用双二次样条模型进行处理。 图 6 拟临界温度示意图 2.2.2 拟临界温度方程的拟合 一定压力对应的拟临界温度是指该压力下等压比热Cp最大时的温度。图 7 为压力分 别为 25 MPa、28 MPa、30 MPa时定压比热Cp随温度的变化曲线。由此图可以看出，在 一定的压力下，随着温度的升高Cp的值开始增加并且在某一温度下达到峰值，然后开始 下降，这个峰值对应的温度即为拟临界温度。在近临界区内，随着压力的增加，拟临界 温度逐渐升高并且Cp对应峰值不断下降。 图 7 不同压力下等压比热随温度变化曲线 在拟临界温度处,有 0=∂ ∂ p Cp τ (4) 本文应用 IAPWS-95 公式以及公式(4)，采用二分法，逐步迭代，精确地计算出了一 系列的拟临界温度值，并拟合了拟临界温度线，方程形式为 6 次多项式： 65432 gxfxexdxcxbxaTcri ++++++= (5) 式中，a、b、c、d、e、f、g为待定系数[8]。对于每一个以拟临界温度线为非规则边界的 计算区域，本文首先确定该区域的压力范围，然后通过迭代计算得到一系列该压力范围 内的拟临界温度值。利用这些值就可以拟合出公式(5)形式的拟临界温度方程，但该方程 只能应用于其对应的非规则区域。 3 结果和讨论 3.1 计算速度 在英特尔酷睿 2 双核 2.4 GHz 处理器和 Microsoft Windows XP 操作系统环境下，通 过 C++语言实现模型的数学计算。对样条插值算法与 IAPWS-IF97 分别在计算时间和计 算时间进行了对比测试，在不同区域均随机了 1000000 个点进行测试。 表 1 ρ(p,T)计算速度（μs） 计算所处区域 样条函数 IAPWS-IF97 规则区域 0.217 IAPWS-IF97：1 区 饱和线附近 0.968 11.343 规则区域 0.217 IAPWS-IF97：2 区 饱和线附近 0.968 13.702 规则区域 0.217 饱和线附近 0.968 IAPWS-IF97：3 区 近临界区 1.079 117.296 图 8 IAPWS-IF97 一区和二区ρ(p,T)与样条函数ρS(p,T)与计算速度比值 图 9 IAPWS-IF97 三区迭代公式ρ(p,T)与样条函数ρS(p,T)与计算速度比值 将发展出的不同样条函数分别与IAPWS-IF97 公式进行计算时间比较，表 1给出了 ρ(p,T)及IAPWS-IF97 各函数的平均计算时间，图 和图 给出了与IAPWS-IF97 的计算 时间比。计算表明样条插值模型ρS(p,T)函数在 1、2 区计算速度是IAPWS-IF97 的 11.7 倍 到 63.1 倍。样条函数在网格相同处理区域均采用统一处理方法，所以在计算不同区域时， 网格处理方法相同的区域计算速度是一样。样条模型规则区域计算速度较快，在非规则 区域要进行网格变形，速度相对降低。在 3 区计算速度是IAPWS-IF97 的 108.7 倍到 540.5 倍，3 区IAPWS-IF97 计算ρS(p,T)时需要迭代，消耗大量时间。 3.2 计算精度 在计算区域内分别测试了样条函数与 IAPWS-95 的最大偏差，和 IAPWS-IF97 与 IAPWS-95 的最大偏差。偏差测试是分为 5 个区进行的，区域的划分如下，采用双精度 计算进行测试。 表 2 样条函数与 IAPWS-95 最大偏差、IAPWS-IF97 与 IAPWS-95 最大偏差 Δ ρ / 最大相对误差 样条函数 IAPWS-IF97 区域 1 4.1457e-6 3.6723e-5 区域 2 8.4863e-6 3.7547e-4 区域 3 7.4682e-6 0.0032148 区域 4 1.9487e-2 0.029573 区域 1： IAPWS-IF97 模型 1 区； 区域 2：除区域 4 包含的所有 IAPWS-IF97 模型 2 区； 区域 3：除区域 4 包含的所有 IAPWS-IF97 模型 3 区； 区域 4：p=20-30 MPa，T=645-680 K。 每个区域随机选取了 1000000 个点进行测试，计算平均偏差如表 2所示。计算结果 表明：以IAPWS-95 为标准数据源，除临界点附近精度与IAPWS-IF97 相当外，在其它区 域偏差均小于IAPWS-IF97，达到很高精度，能够满足工业生产的要求。 4 结论 a) 本文引入双二次样条插值模型应用于水和水蒸气密度计算，计算速度是 IAPWS-IF97 的 11.7 倍到 540.5 倍。在 IAPWS-IF97 模型 1、2 区，IAPWS-IF97 以 p， T 作为自变量的比吉布斯自由能 g(p,T)作为基本函数，计算速度较快，样条插值模 型计算速度是它的 11.7 倍到 63.1 倍。在 3 区，IAPWS-IF97 以亥姆霍兹自由能 f(ρ,T) 为基本函数，以 p，T 为自变量求解其它热物性参数时需要迭代，样条模型的计算 速度是它的 108.7 倍到 540.5 倍。 b) 样条插值模型以 IAPWS-95 为数据源进行插值，除临界点附近外，计算精度均要高 于 IAPWS-IF97，能够满足工业生产的实际需求。 c) 相对于 IAPWS-IF97 的分区，样条插值模型可以用于跨区域计算；避免了分区带来 的区域边界不一致性问题，并且能够连续的表征热力学面。 d) 样条插值模型实现了以 p，T 为自变量密度函数的快速精确计算，可以进一步应用 于其他热物性为自变量的计算。 参考文献 [1] Wagner W and Cooper J R, The IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam. Journal of Engineering for Gas Turbines and Power, ASME, 2000, 122: 150-182 [2] Miyagawa, K and Hill, P G., A Tabular Taylor Series Expansion Method for Fast Calculation of Steam Properties. Journal of Engineering for Gas Turbines and Power, ASME, 1997, 119: 485-491 [3] Miyagawa K and Hill P G , Rapid and Accurate Calculation of Water and Steam Properties Using the Tabular Taylor Series Expansion Method, Journal of Engineering for Gas Turbines and Power, ASME, 2001, 123: 707-712 [4] Späth H, Two-Dimensional Spline Interpolation Algorithms, USA: Wellesley, MA, 1995: 31– 69 [5] Kunick M, Kretzschmar, H J and Gampe U , Fast Calculation of Thermodynamic Properties of Water and Steam in Process Modelling using Spline Interpolation, 15th International Conference on the Properties of Water and Steam, Berlin. 2008 [6] Wagner W and Pruß A, The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use, J. Phys. Chem. Ref. Data, 2002,31: 387 – 485 [7] 晏永飞，陈保东，王伟. 基于IAPWS-IF97的第3区水和水蒸汽密度的迭代计算.辽宁石油化工大学 学报, 2007, 27(4): 38 – 40 Yan Y F, Chen B D and Wang W, Iterative Computation of Water and Steam Density Based on t he Third Zone of IAPWS-IF97, Journal of Liaoning University of Petroleum & Chemical Technology, 2007, 27(4): 38 – 40 [8] 梁法春，王栋，林宗虎. 超临界区水的你临界温度的确定. 动力工程, 2004, 24(6): 869 – 870. [9] Liang F C, Wang D and Lin Z H, Determination of Pseudocritical Temperture of Water in the Supercritical Pressure Region”, Power Engineering, 2004, 24(6): 869 – 870 水和水蒸气密度的跨区域快速计算方法 作者： 安宾， 王智学， 王晓东， 段远源 作者单位： 安宾,王智学,王晓东(北京科技大学机械工程学院热能工程系,北京,100083)， 段远源(清华大学热科学与动 力工程教育部重点实验室,北京,100084) 本文链接：http://d.g.wanfangdata.com.cn/Conference_7469881.aspx