Fluent表达式及UDF求解温度反演问题

详细信息

作者优秀: 优秀教师/意见领袖/博士学历/特邀专家
平台推荐: 内容稀缺

4月前浏览6406

本案例演示在Fluent中利用表达式功能及UDF处理温度反问题的基本思路。

1 问题描述

案例要解决的问题如下图所示。管道直径20 mm，长度300 mm。流体介质为空气，入口速度2 m/s，出口静压0 Pa，固体壁面对流换热系数 20 W/(m2-K)，环境温度300 K，若想要管道出口平均温度为350 K，请问管道入口温度应为多少？

本案例几何简单，可以利用传热学的相关理论直接计算。不过这仅仅只是为了方便演示而已，如果计算模型非常复杂的话，利用传热学公式也不太容易计算。

这是个典型的反问题，利用仿真的方式也不容易处理。一种解决方式是采用二分法不断尝试入口温度，使得出口平均温度接近目标值。如本案例中，可以先给入口一个温度上限（比如上限600 K），然后用中间值500 K及上限值600 K分别进行计算并查看出口值，若600 K入口计算得到的出口温度值大于350 K，而500 K温度值计算得到的温度值也大于350 K，则使用350~500的中间值（425 K）进行计算并比较出口温度值；若500 K的温度值计算得到的出口温度值小于350 K，则使用500~600的中间值（550 K）进行计算并比较出口值。相同的方式不断逼近，得到最终的入口温度值。另一种方式也可以不采用二分法，而是采用单步搜索方式对入口温度进行修改。

两种方法在技术上都是可以行得通的，然而在具体实施上不太方便，需要多次修改边界条件并执行计算以及查看出口温度并获取入口温度。

本案例借助Fluent中的表达式功能将出口温度作为入口温度的操作变量，在迭代计算过程中自动修改入口温度，使出口温度满足设计目标。

注：本案例也可以使用UDF中的DEFINE_ADJUST及DEFINE_PRFILE宏进行处理。
”

2 Fluent设置

进入Fluent并导入网格

2.1 Models设置

激活能量方程

使用默认的SST k-omega湍流模型

2.2 Materials设置

材料介质air采用默认参数

2.3 设置边界条件

1、inlet边界

设置入口速度为2 m/s

设置入口温度350 K

注：这里的入口温度是随便给的，后面要重新指定。
”

2、outlet边界

指定出口静压0 Pa

指定出口回流温度300 K

3、wall边界

指定壁面换热系数20 W/(m2-K)，指定环境温度300 K

2.4 定义报告

如下图所示添加报告

报告名称为report-outlet-temp，指定其为出口平均温度

2.5 初始化

进行初始化

2.6 进行计算

设置迭代300步进行计算

监测得到的出口平均温度为336.82 K，而需要的出口目标温度为350 K，显然入口温度需要调整。

3 表达式处理

3.1 定义表达式

右键选择模型树节点Named Expressions，点击弹出菜单项**New…**新建表达式

如下图所示创建表达式Tout，其定义为前面所指定的出口温度报告值report-outlet-temp

定义表达式TInCurrent，其为入口边界的当前值，如下图所示

定义表达式relax，指定其值为0.2，如下图所示

注：这是一个松弛因子，会影响到计算的收敛性。如果后面的计算出现较大的振荡，建议降低此参数值。
”

定义表达式TAdjust，如下图所示

这里350 K为目标温度，这里在入口当前温度值的前后尝试
”

定义表达式TIn，如下图所示

3.2 修改边界条件

指定入口边界inlet的温度为TIn，如下图所示

3.3 定义报告

如下图所示，新建一个Expression报告

指定其值为TIn，如下图所示

3.4 进行计算

指定初始温度为300 K

注：这里将初始温度设置为300 K仅仅是为了方便观察出口温度分布。实际上设置初始温度为350 K要更好一些，计算收敛更快。
”

设置迭代计算300步

3.5 计算结果

出口温度的变化曲线

查看入口温度监测曲线，可以看到最终温度为367.35 K，即入口温度在367.35 K时，出口温度能稳定在350 K。

4 UDF处理

案例所涉及的问题也可以使用UDF进行处理，其处理思路与表达式相同。

4.1 编译并加载UDF

编写UDF源代码如下。

#include "udf.h"

#define ZoneID 6      //定义为出口的ID，后面统计出口温度时要用
#define Target_T 350  //定义你希望达到的出口温度
#define Initial_T 350 //定义入口温度的初始值

// 定义全局变量，便于宏内调用
real staT = 0.0;           //初始化为0
real currentT = Initial_T; //初始化为温度初始值

// 计算入口的温度值
DEFINE_ADJUST(in_avg_T, domain)
{
    real NV_VEC(A);
    real sumTA = 0.0, sumA = 0.0; //声明实数类型变量

#if !RP_HOST
    face_t f;
    Thread *fthread;
    fthread = Lookup_Thread(domain, ZoneID);
    begin_f_loop(f, fthread)
    {
        if (PRINCIPAL_FACE_P(f, fthread))
        {
            F_AREA(A, f, fthread);
            sumA  = NV_MAG(A);
            sumTA  = NV_MAG(A) * F_T(f, fthread);
        }
    }
    end_f_loop(f, fthread)
#endif
        // 变量全局约简并从计算节点传递到host节点
        sumTA = PRF_GRSUM1(sumTA);
    sumA = PRF_GRSUM1(sumA);
    node_to_host_real_2(sumTA, sumA);

#if RP_HOST
    //得到温度的面平均值
    staT = sumTA / sumA;
    Message("the Temperature at outlet is %e K. \n ", staT);
#endif
    host_to_node_real_1(staT);
}

// 实现入口温度在计算过程中逐渐调节，采用比例调节，比例系数0.2
// 比例系数取值范围0~1，调节方式等同于亚松弛因子
DEFINE_PROFILE(inlet_temperature, thread, i)
{
    face_t f;
    //迭代次数大于5 时且每隔10次迭代更新一次入口温度
    if (N_ITER > 5 && N_ITER % 10 == 0)
    {
        //当出口温度值与目标温度值差异超过0.1K时，更新入口温度
        if (abs(staT - Target_T) > 0.1)
        {
            currentT = currentT - (0.2 * (staT - Target_T));
            //将计算出的入口温度值赋值在入口边界上
            begin_f_loop(f, thread)
            {
                F_PROFILE(f, thread, i) = currentT;
            }
            end_f_loop(f, thread)
        }
    }
    else
    {
        begin_f_loop(f, thread)
        {
            F_PROFILE(f, thread, i) = currentT;
        }
        end_f_loop(f, thread)
    }
}