Fluent中调用IAPWS-IF97材料模型

在CFD计算过程中，材料属性参数是直接影响计算结果的重要输入数据之一。对于水/水蒸气介质参数，目前工程上普遍使用IAPWS-IF97模型。本文描述在Fluent中利用UDF调用IAPWS-IF97材料数据模型。

Fluent中并没有集成IAPWS-IF97，只提供了NIST以及真实气体模型，然而如果是处理水与水蒸气，还是IAPWS-IF97比较好用。Fluent中的NIST模型毛病太多，而真实气体模型又太过于复杂。虽然利用Fluent提供的UDRGM通过编程可以实现IAPWS-IF97的效果，但对于大多数人来讲实在是强人所难了。

提示：CFX中集成了IAPWS-IF97，如果不想折腾的话，不妨使用CFX。

IAPWS-IF97有python版本接口（https://github.com/jjgomera/iapws），不过如果能有C/C++版本的话，利用UDF调用起来也相对简单一点。如果用Python，需要先打包成dll后通过UDF调用，中间过程极其复杂且容易出错，想想还是算了，本来利用商业软件就是为了省事儿。

下面是本文的主角：IF97。官网地址：https://github.com/CoolProp/IF97，出自之前的文章介绍过的CoolProp（流体热物性计算库CoolProp）。其为原生C++语言编写，所有内容只有一个名为IF97.h的头文件，调用的时候直接包含即可，非常的简单。关于IF97的使用，官网上提供了一个名为IF97.cpp的文件，其中包含了函数调用过程。当然自己分析代码也很容易找到函数的使用方式。

使用的时候注意温度和压力范围：

或

注：更高的压力和温度怎么办？此时IAPWS-IF97失效了，没有更好的办法了，试验数据拟合吧。

使用的时候注意单位，默认情况下使用国际单位制。若在源文件中指定了#define IAPWS_UNITS，则压力使用MPa，所有能量单位为kJ。

一个简单的示例如下：

#include "IF97.h"
#include 
using namespace std;
int main() {
    using namespace IF97;
    //温度500K，压力3MPa
    double T = 500, p = 3e6;
    cout << "-------------T = 500K, p = 3e6 Pa----------------------" << endl;
    cout << "密度:" << rhomass_Tp(T, p) << "kg / m3" << endl;
    cout << "等压比热容:" << cpmass_Tp(T, p) << " J / (kg - K)"<< endl;
    cout << "等容比热容:" << cvmass_Tp(T, p) << " J / (kg - K)"<< endl;
    cout << "声速:" << speed_sound_Tp(T, p)<<"m/s" << endl;
    cout << "质量焓:" << hmass_Tp(T, p) << "J / kg"<<endl;
    cout << "质量熵:" << smass_Tp(T, p) << "J / (kg - K"<< endl;
    cout << "质量内能:" << umass_Tp(T, p) << "J / kg"<< endl;
    cout << "动力粘度:" << visc_Tp(T, p) << "Pa.s" << endl;
    cout << "热导率:" << tcond_Tp(T, p) << "W/(m-K)" << endl;
    cout << "prandtl数:" << prandtl_Tp(T, p) << endl;
    cout << "-------------------------------------" << endl;
 
 
    //温度1500K，压力30MPa
    double T1 = 1500, p1 = 30e6;
    cout << "-------------T =1500K, p = 30e6 Pa----------------------" << endl;
    cout << "密度:" << rhomass_Tp(T1, p1) << "kg / m3" << endl;
    cout << "等压比热容:" << cpmass_Tp(T1, p1) << " J / (kg - K)" << endl;
    cout << "等容比热容:" << cvmass_Tp(T1, p1) << " J / (kg - K)" << endl;
    cout << "声速:" << speed_sound_Tp(T1, p1) << "m/s" << endl;
    cout << "质量焓:" << hmass_Tp(T1, p1) << "J / kg" << endl;
    cout << "质量熵:" << smass_Tp(T1, p1) << "J / (kg - K" << endl;
    cout << "质量内能:" << umass_Tp(T1, p1) << "J / kg" << endl;
    cout << "动力粘度:" << visc_Tp(T1, p1) << "Pa.s" << endl;
    cout << "热导率:" << tcond_Tp(T1, p1) << "W/(m-K)" << endl;
    cout << "prandtl数:" << prandtl_Tp(T1, p1) << endl;
    cout << "-------------------------------------" << endl;
    return 0;
}

程序输出为：

------T = 500K, p = 3e6 Pa------------
密度:831.658kg / m3
等压比热容:4655.81 J / (kg - K)
等容比热容:3221.39 J / (kg - K)
声速:1240.71m/s
质量焓:975542J/kg
质量熵:2580.42J/ kg-K)
质量内能:971935J/kg
动力粘度:0.000117996Pa.s
热导率:0.63979W/(m-K)
prandtl数:0.858669
-------------------------------------
--------T =1500K, p = 30e6 Pa--------
密度:43.3348kg / m3
等压比热容:2727.24 J/(kg-K)
等容比热容:2192.75 J/(kg-K)
声速:928.548m/s
质量焓:5.16724e+06J/kg
质量熵:7729.7J/(kg-K)
质量内能:4.47495e+06J/kg
动力粘度:5.69793e-05Pa.s
热导率:0.198873W/(m-K)
prandtl数:0.781384
-------------------------------------

如下图所示。

一些常见的热物性及对应的函数如下表所示。

在UDF中调用IF97也非常方便，只不过需要注意的是，IF97是采用C++语言编写，编译的时候需要注意。

下面随便编写一个UDF进行测试。

#include "udf.h"
#include "IF97.h"
using namespace IF97;
DEFINE_ON_DEMAND(testIF97)
{
    real T= 500;
    real p = 3e6;
    printf("rho:%gn",rhomass_Tp(T,p));
}

这里利用GCC编译工具进行编译。可以采用之前公众号推出的利用GCC进行编译的方式，这里图省事儿，利用硫酸亚铜提供的编译工具（见网址http://blog.sina.com.cn/s/blog_14d64daa10102z7wy.html），编译之前记得将源文件和头文件的编码格式都修改为ANSI（默认为UTF-8），否则编译会出现各种莫名其妙的错误信息。GCC就是矫情！

• 编译完毕后加载UDF

• 执行UDF

• UDF运行后的结果如下图所示

可以看到，UDF正确的计算出了温度500 K，压力3 MPa条件下的水的密度、比热、粘度以及热导率。

好了，扯了那么多，自然不是为了在屏幕上打印一堆没用的信息，我们需要利用UDF宏调用IF7计算材料介质物理属性，这里需要编写DEFINE_PROPERTY宏。CoolProp提供了一些调用的UDF宏，但是只提供了Linux下的编译方式，如何在Windows下调用，目前还有待摸索。

#include "udf.h"
#include "IF97.h"
using namespace IF97;
using namespace std;
const real gauge = 101325;/*操作压力*/
 
//定义密度
DEFINE_PROPERTY(water_density,c,t)
{
    real temperature = C_T(c,t);
    real pressure = C_P(c,t) + gauge;
    real density = rhomass_Tp(temperature,pressure);
    return density;
}
 
//定义粘度
DEFINE_PROPERTY(water_viscosity,c,t)
{
    real viscosity;
    real temperature = C_T(c,t);
    real pressure = C_P(c,t) + gauge;
    real density = rhomass_Tp(temperature,pressure);
    viscosity = visc_Tp(temperature,pressure);
    return viscosity;
}
 
//定义热导率
DEFINE_PROPERTY(water_thermalConductivity,c,t)
{
    real thermalConductivity;
    real temperature = C_T(c,t);
    real pressure = C_P(c,t) + gauge;
    real density = rhomass_Tp(temperature,pressure);
    thermalConductivity = tcond_Tp(temperature,pressure);
}
 
//定义比热容
DEFINE_SPECIFIC_HEAT(water_specificHeat,temperature,Tref,enthalpy,yi)
{
    real pressure;
    real specificHeat;
    Domain *domain = Get_Domain(1);
    Thread *t;
    cell_t c;
 
    thread_loop_c(t,domain)
    {
        begin_c_loop(c,t)
        {
            pressure = C_P(c,t)+gauge;
            temperature = C_T(c,t);
            specificHeat = cpmass_Tp(temperature,pressure);
            *enthalpy = specificHeat*(temperature-Tref);
            return specificHeat;
        }end_c_loop(c,t)
    }
}
 
//这段代码没有用
DEFINE_ON_DEMAND(testIF97)
{
    real T= 500;
    real p = 3e6;
 
    #if RP_HOST
    Message("water Properties(T=500K,p=3MPa)n");
    Message("-----------------------------n");
    Message("rho:%gn",rhomass_Tp(T,p));
    Message("cp:%gn",cpmass_Tp(T,p));
    Message("viscosity:%gn",visc_Tp(T,p));
    Message("Thermal Conductivity:%gn",tcond_Tp(T,p));
    #endif
}

• 编译后加载UDF

• 加载后控制台窗口提示如下图所示

• 此时可以在材料对话框中加载UDF了

随便弄个案例测试一下。如下弯管，入口速度10 m/s，温度400 K，出口静压0 Pa，所有壁面绝热，计算管内流场。

• 压力分布

• 密度分布

• 速度分布

建议各位道友自己尝试着编译UDF，若实在搞不定的话，花点银子直接下载下面的文件吧，毕竟有钱能使鬼推磨。

案例相关文件下载链接：

本篇文章来源于微信公众号: CFD之道