前面的文章中介绍了利用IF97开源代码^[1]进行水-水蒸气材料介质属性定义。然而此代码是利用C++语言编写，常规情况下编译不太容易，本文将其改写为C代码形式。

IAPWS-IF97^[2]中饱和压力与饱和温度的计算位于第4区。

描述饱和线的方程是一个隐式二次方程，该方程可以直接用于求解饱和压力  与饱和温度 。该方程描述为：

式中，

其中，，，参数  如下表所示。

1 饱和压力方程

通过式（1）可求解出饱和压力计算公式为：

式中，，且有：

式中  为方程（3）计算得到。

饱和压力方程的适用范围：

2 饱和温度方程

饱和温度计算公式为：

式中，，且有：

其中：

参数  及  在前面已经给出。饱和温度方程的使用范围：

整理成代码：

/*IF97.H*/
#include "math.h"
#include "stdio.h"
 
/*设置n[0]=0仅为了后面公式中ni编号一致*/
static double n[] = {
    0,
    0.11670521452767e4,
    -0.72421316703206e6,
    -0.17073846940092e2,
    0.12020824702470e5,
    -0.32325550322333e7,
    0.14915108613530e2,
    -0.48232657361591e4,
    0.40511340542057e6,
    -0.23855557567849,
    0.65017534844798e3,
};
/*根据温度计算饱和压力，温度单位为K，压力单位Pa*/
double p_T(double T)
{
    double T_star = 1;   /*默认为1 K*/
    double p_star = 1e6; /*默认为1 MPa*/
    if ((T > 647.096) || (T < 273.15))
    {
        printf("Temperature out of range");
        return 0.0;
    }
    else
    {
        double theta = T / T_star + n[9] / (T / T_star - n[10]);
        double A = theta * theta + n[1] * theta + n[2];
        double B = n[3] * theta * theta + n[4] * theta + n[5];
        double C = n[6] * theta * theta + n[7] * theta + n[8];
        return p_star * pow(2 * C / (-B + sqrt(B * B - 4 * A * C)), 4);
    }
}
/*根据压力计算饱和文件，压力单位Pa，温度单位K*/
double T_p(double p)
{
    double p_star = 1e6;
    double T_star = 1;
 
    if ((p < 611.213) || (p > 22.064E6))
    {
        printf("Pressure out of range");
        return 0.0;
    }
    else
    {
        double beta2 = sqrt(p / p_star);
        double beta = sqrt(beta2);
        double EFG[3];
        double &E = EFG[0];
        double &F = EFG[1];
        double &G = EFG[2];
 
        EFG[0] = 1.0;
        EFG[1] = n[1];
        EFG[2] = n[2];
        for (int i = 0; i < 3; ++i)
        {
            EFG[i] *= beta;
        }
 
        for (int i = 0; i < 3; ++i)
        {
            EFG[i] += n[i + 3];
        }
 
        for (int i = 0; i < 3; ++i)
        {
            EFG[i] *= beta;
        }
 
        for (int i = 0; i < 3; ++i)
        {
            EFG[i] += n[i + 6];
        }
 
        const double D = 2 * G / (-F - sqrt(F * F - 4 * E * G));
        const double n10pD = n[10] + D;
        return T_star * 0.5 * (n10pD - sqrt(n10pD * n10pD - 4 * (n[9] + n[10] * D)));
    }
}

随便写个文件测试一下。

/*IF97.CPP*/
#include "IF97.h"
#include 
using namespace std;
 
int main()
{
    double T[3] = {300, 500, 600};
    double p[3] = {1e5, 1e6, 10e6};
    for (int i = 0; i < 3; i++)
    {
        cout << "Temperature:" << T[i] << "K,Psat:" << p_T(T[i]) << " Pa" << endl;
    }
    for (int i = 0; i < 3; i++)
    {
        cout << "Pressure:" << p[i] << " Pa,Tsat:" << T_p(p[i]) << " K" << endl;
    }
    return 0;
}

输出结果：

Temperature:300K,Psat:3536.59 Pa
Temperature:500K,Psat:2.6389e+06 Pa
Temperature:600K,Psat:1.23443e+07 Pa
Pressure:100000 Pa,Tsat:372.756 K
Pressure:1e+06 Pa,Tsat:453.036 K
Pressure:1e+07 Pa,Tsat:584.149 K

前面的代码可以直接拷贝到UDF中使用。

历史文章：Fluent中调用IAPWS-IF97材料模型

参考资料