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

FLUENT UDF 应用实例:传热问题第二第三类热边界条件转
换成第一类边界条件
1 引言
传热问题的常见边界条件可归纳为三类，以稳态传热为例，三类边界条件的表达式如下。

恒温边界（第一类边界条件）：const w T = （1-1）
恒热流密度边界（第二类边界条件）：const w T n λ∂⎛⎫-= ⎪∂⎝⎭ （1-2）
对流换热边界（第三类边界条件）：()w f w
T h T T n λ∂⎛⎫
-=- ⎪∂⎝⎭ （1-3）
2 问题分析
2.1 纯导热问题
以二维稳态无源纯导热问题为例，如图1所示，一个10×10m 2的方形平面空间，上下面以及左边为恒温壁面（21℃），右边第二类、第三类边界条件如图所示。

为方便问题分析，内部介质的导热系数取1W/m ℃。

模型水平垂直方向各划分40个网格单元，不计边界条件处壁厚。

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

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

第二类边界条件：11M M q T T δ
λ
-=+
（2-1）
第三类边界条件：11/1M M f
h h T T T δδ
λλ-⎛
⎫⎛⎫
=++
⎪ ⎪⎝
⎭⎝
⎭
（2-2） 其中，T M 为边界节点处的温度（所求值），T M-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)
}。