膜态沸腾UDF【转载】

膜态沸腾的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)

    }

}