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代码
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); //获取边界邻近网格单元的IDF_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); //获取边界邻近网格单元的IDF_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)}
