FLUENT UDF应用实例:传热问题第二第三类热边界条件转换成第一类边界条件

1 引言

    传热问题的常见边界条件可归纳为三类，以稳态传热为例，三类边界条件的表达式如下。

    恒温边界（第一类边界条件）：                       （1-1）

    恒热流密度边界（第二类边界条件）：   （ 1-2）

    对流换热边界（第三类边界条件）：   （1-3）

2 问题分析

2.1 纯导热问题

    以二维稳态无源纯导热问题为例，如图1所示，一个10×10m2的方形平面空间，上下面以及左边为恒温壁面（21℃），右边第二类、第三类边界条件如图所示。为方便问题分析，内部介质的导热系数取1W/m℃。模型水平垂直方向各划分40个网格单元，不计边界条件处壁厚。

图1 问题描述

采用FLUENT软件自带边界条件直接进行计算，结果如图2所示。

（a）第二类边界条件

（b）第三类边界条件

图2 软件自带边界计算结果

参考数值传热学[3]，对于第二类（式1-2）、第三类（式1-3）边界条件可通过补充边界点代数方程的方法进行处理，结果如下。

第二类边界条件：                （2-1）

第三类边界条件：     （2-2）

    其中，TM为边界节点处的温度（所求值），TM-1为靠近边界第一层网格节点处的温度， 为靠近边界第一层网格节点至边界的法向距离，q为热流密度，h为对流换热系数。

将以上两式通过UDF编写成边界条件（DEFINE_PROFILE），全部转换为第一类边界条件，计算结果如图3所示。

（a）第二类边界条件

（b）第三类边界条件

图3 UDF计算结果

可以看出，经过UDF边界转换后的计算结果与软件自带边界计算结果几乎完全相同。

2.2 对流换热问题

    以上处理方式对于导热问题肯定是适用的，但是对于对流换热问题能否用同样的方式处理呢，笔者认为，严格意义上讲式2-1和2-2对与对流换热问题是不能用的，因为边界内侧的流体与壁面的换热机制是对流换热。但是，当第一层网格重心距离边界很近时，该位置流体的流速很小，无限接近边界时流速应该为零。因此只要第一层网格重心足够靠近边界，换热机制用导热来处理应该是没有问题的。为此，建立图4所示的计算模型，尺寸、网格、介质同上，此时首层网格重心与边界的法向距离为0.125m。

图4 问题描述

采用FLUENT软件自带边界条件直接进行计算，结果如图5所示。

（a）第二类边界条件

（b）第三类边界条件

图5 软件自带边界计算结果

采用UDF将第二类、第三类边界条件转换成第一类边界条件计算结果如图6所示。

（a）第二类边界条件

（b）第三类边界条件

图6 UDF计算结果

对比可知，首层网格重心与边界的法向距离为0.125m时，UDF计算结果与软件自带边界条件计算结果差距较大。因此，对壁面处网格进行细化如图7所示，此时首层网格重心与边界的法向距离为0.0625m，重新计算，结果如图8、图9所示。可以看出，细化网格后，经过UDF边界转换后的计算结果与软件自带边界计算结果几乎完全相同。说明在该厚度内，导热成为传热的主导机制，对流换热几乎可以忽略了，对于不同的问题，这个厚度到底达到多少也不确定，但是可以确定的是至少需要在边界层内，越靠近壁面，这种处理越准确。继续进行细化，两者的计算结果同样几乎一致。

图7 壁面网格细化

（a）第二类边界条件

（b）第三类边界条件

图8 软件自带边界条件计算结果

（a）第二类边界条件

（b）第三类边界条件

图9 UDF计算结果

3 总结

    换热问题的第二类、第三类边界条件都可以同UDF转换成第一类边界条件进行计算，但是边界处网格需进行细化，特别对于对流换热问题，首层网格重心与边界的法向距离至少要在边界层范围内，距离越小，转换计算的准确性越高。

参考文献：

[1]FLUENT帮助文档

[2]杨世铭,陶文铨.传热学(第4版)[M].北京:高等教育出版社

[3]陶文铨.数值传热学[M].西安:西安交通大学出版社

附录 UDF代码

#include "udf.h"
#include "prop.h"
#include "mem.h"
DEFINE_SOURCE(uds_source,c,t,dS,eqn)
{
double source=10;
dS[eqn]=0;
return source;
}
DEFINE_PROFILE(MyProfile2,thread,index)
{
double q=10.0;
double T;
face_t f;
cell_t c0;
Thread *t0;
cell_t c; 
double thickness;
double thick;
double x[ND_ND];
double x0[ND_ND];
double NV_VEC(a);
double NV_VEC(A);
double cosphi; 
begin_f_loop(f,thread) 
{
t0=THREAD_T0(thread); //获取边界邻近网格的线
c0=F_C0(f,thread);  //获取边界邻近网格单元的ID
F_CENTROID(x,f,thread); //获取面心
C_CENTROID(x0,c0,t0);   //获取邻近单元重心
NV_VV(a, =, x0, -, x); //求面心和重心的向量(两心向量）
thick=NV_MAG(a);  //求面心和重心的距离（非法向）
F_AREA(A,f,thread);  //获取边界法向面积矢量
cosphi=fabs(NV_DOT(a,A))/NV_MAG(a)/NV_MAG(A); //求两心向量和法向的夹角余弦
thickness=thick*cosphi;  //求两心法向距离
T=C_T(c0,t0)+q*thickness;
F_PROFILE(f,thread,index)=T;    
//printf("\n 首层网格重心到边界的法向距离为 %g\n",thickness);
}   
end_f_loop(f,thread)
}
DEFINE_PROFILE(MyProfile3,thread,index)
{
double coef=5.0;
double tout=10;
double T;
face_t f;
cell_t c0;
Thread *t0;
cell_t c; 
double thickness;
double thick;
double x[ND_ND];
double x0[ND_ND];
double NV_VEC(a);
double NV_VEC(A);
double cosphi; 
begin_f_loop(f,thread) 
{
t0=THREAD_T0(thread); //获取边界邻近网格的线
c0=F_C0(f,thread);  //获取边界邻近网格单元的ID
F_CENTROID(x,f,thread); //获取面心
C_CENTROID(x0,c0,t0);   //获取邻近单元重心
NV_VV(a, =, x0, -, x); //求面心和重心的向量(两心向量）
thick=NV_MAG(a);  //求面心和重心的距离（非法向）
F_AREA(A,f,thread);  //获取边界法向面积矢量
cosphi=fabs(NV_DOT(a,A))/NV_MAG(a)/NV_MAG(A); //求两心向量和法向的夹角余弦
thickness=thick*cosphi;  //求两心法向距离
T=(C_T(c0,t0)+coef*thickness*(tout+273.15))/(1+coef*thickness);
F_PROFILE(f,thread,index)=T;    
//printf("\n 首层网格重心到边界的法向距离为 %g\n",thickness);
}   
end_f_loop(f,thread)
}