膜态沸腾的UDF,添加注释。其中获取VOF梯度的方法详见前面的日志,其中很多宏无法通过UDF手册查阅,
蒸汽相中的质量源项的一般形式为:
式中:
通过一阶近似,热流之差可表达为:
式中:
通过此近似,源项变为:
因为没有内部质量源项,所以也想的质量源项变为:
能源方程的潜热变为:
截面的物性包括,表面张力为0.1N/m,潜热为105J/kg,饱和温度为500K。该问题的长度为Taylor-Raleigh不稳定性的波长。
该问题的速度为:
因此时间尺度为:
计算域水平宽度为,垂直高度为
。网格尺度为64(水平方向)×192(垂直方向),蒸汽液体截面的初始形状受到气泡增长的影响。因此,需要另外一个初始化UDF
式中,x(y)为水平(垂直)轴,单位为m。
努塞尔数为表征沸腾换热的无量纲数,其定义为:
因为该问题的时间为0.01s,时间步长为2×10-4,50个时间步。该问题需要进行6000步时间步迭代来捕捉气泡释放的过程
源代码如下:
#include "udf.h"
#include "sg.h"
#include "sg_mphase.h"
#include "flow.h"
#include "mem.h"
DEFINE_ADJUST(area_density, domain)
{
Thread *t; //定义一个线程指针
Thread **pt; //定义一个指向线程指针的指针
cell_t c; //定义单元
Domain *pDomain = DOMAIN_SUB_DOMAIN(domain,P_PHASE);
real voidx, voidy, voidz=0;
{
Alloc_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
Scalar_Reconstruction(pDomain, SV_VOF,-1,SV_VOF_RG,NULL);
Scalar_Derivatives(pDomain,SV_VOF,-1,SV_VOF_G,SV_VOF_RG,
Vof_Deriv_Accumulate);
}
{
Alloc_Storage_Vars(domain, SV_T_RG, SV_T_G, SV_NULL);
T_derivatives(domain);
Free_Storage_Vars(domain, SV_T_RG, SV_NULL);
}
mp_thread_loop_c (t,domain,pt)
if (FLUID_THREAD_P(t))
{
Thread *tp = pt[P_PHASE];
begin_c_loop (c,t)
{
#if RP_3D
C_UDMI(c,t,0) = (C_VOF_G(c,tp)[0]*C_T_G(c,t)[0]+
C_VOF_G(c,tp)[1]*C_T_G(c,t)[1]+C_VOF_G(c,tp)[2]*C_T_G(c,t)[2]);
#endif
#if RP_2D
C_UDMI(c,t,0) = (C_VOF_G(c,tp)[0]*C_T_G(c,t)[0]+
C_VOF_G(c,tp)[1]*C_T_G(c,t)[1]);
#endif
}
end_c_loop (c,t)
}
Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
Free_Storage_Vars(domain, SV_T_G, SV_NULL);
}
DEFINE_SOURCE(gas, cell, thread, dS, eqn)
{
real x[ND_ND]; //定义一个变量用来存储网格的位置信息
real source;
Thread *tm = THREAD_SUPER_THREAD(thread); //获取混合相的指针
Thread **pt = THREAD_SUB_THREADS(tm);
real Kl = C_K_L(cell, pt[1])*C_VOF(cell, pt[1]),
Kg = C_K_L(cell, pt[0])*C_VOF(cell, pt[0]);
real L = 1e5; //水的气化潜热
source = (Kl+Kg)*C_UDMI(cell,tm,0) / L;
C_UDMI(cell, tm, 1) = source;
C_UDMI(cell, tm, 2) = -source*L;
dS[eqn] =0;
return source;
}
DEFINE_SOURCE(liquid, cell, thread, dS, eqn)
{
real x[ND_ND];
real source;
Thread *tm = THREAD_SUPER_THREAD(thread);
Thread **pt = THREAD_SUB_THREADS(tm);
source = -C_UDMI(cell, tm, 1);
dS[eqn] = 0;
return source;
}
DEFINE_SOURCE(energy, cell, thread, dS, eqn)
{
real x[ND_ND];
real source;
Thread *tm = thread;
source = C_UDMI(cell, tm, 2);
dS[eqn] = 0;
return source;
}
DEFINE_INIT(my_init_function, domain)
{
Thread *t;
Thread **pt;
Thread **st;
cell_t c;
Domain *pDomain = DOMAIN_SUB_DOMAIN(domain,P_PHASE); //获取主相的指针
Domain *sDomain = DOMAIN_SUB_DOMAIN(domain,S_PHASE); //获取次相的指针
real xc[ND_ND], y, x;
mp_thread_loop_c (t,domain,pt)
if (FLUID_THREAD_P(t))
{
Thread *tp = pt[P_PHASE]; //获取主相的指针
begin_c_loop (c,t)
{
C_CENTROID(xc,c,t);
x=xc[0];
y=xc[1];
if ( y < 0.00292 + 0.0006*cos(6.283*x/0.0778) )
C_VOF(c,tp) = 1;
else
C_VOF(c,tp) = 0;
}
end_c_loop (c,t)
}
mp_thread_loop_c (t,domain,st)
if (FLUID_THREAD_P(t))
{
Thread *sp = st[S_PHASE]; //获取次相的指针
begin_c_loop (c,t)
{
C_CENTROID(xc,c,t);
x=xc[0];
y=xc[1];
if ( y < 0.00292 + 0.0006*cos(6.283*x/0.0778) )
C_VOF(c,sp) = 0;
else
C_VOF(c,sp) = 1;
}
end_c_loop (c,t)
}
}