膜态沸腾的UDF，添加注释。其中获取VOF梯度的方法详见前面的日志，其中很多宏无法通过UDF手册查阅，

蒸汽相中的质量源项的一般形式为：

式中：

通过一阶近似，热流之差可表达为：

式中：

通过此近似，源项变为：

因为没有内部质量源项，所以也想的质量源项变为：

能源方程的潜热变为：

截面的物性包括，表面张力为0.1N/m，潜热为10⁵J/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)

}

}