日志
2018-04-11
|||
#include "udf.h"
#include "sg_mphase.h"
#include "sg_vof.h"
#include "sg.h"
#include "mem.h"
#include "flow.h"
#include "metric.h"
#define VISCOSITY 0.001
#define SURF_TENS 0.0728
#define MYTRUE 1
#define MYFALSE 0
#define Hoff(x) acos(1-(2.0*tanh(5.16*(pow((x/(1+(1.31*pow(x,0.99)))),0.0706)))))
#define static_Con_Ang 95
#define domain_ID 2
DEFINE_ADJUST(adjust_VOF_gradient,domain)
{
Thread *t;
cell_t c;
face_t f;
domain = Get_Domain(domain_ID);
thread_loop_c(t,domain)
{
begin_c_loop(c,t)
{
C_UDSI(c,t,0)=C_VOF(c,t);
}
end_c_loop(c,t)
}
thread_loop_f(t,domain)
{
if(THREAD_STORAGE(t,SV_UDS_I(0))!=NULL)
begin_f_loop(f,t)
{
F_UDSI(f,t,0)=F_VOF(f,t);
}
end_f_loop(f,t)
}
}
DEFINE_ON_DEMAND(store_gradient)
{
Domain *domain;
cell_t c;
Thread *t;
domain = Get_Domain(1);
thread_loop_c(t,domain)
{
begin_c_loop(c,t)
{
C_UDMI(c,t,0)=NV_MAG(C_UDSI_G(c,t,0));
}
end_c_loop(c,t)
}
}
DEFINE_PROFILE(con_ang,thread,position)
{
Thread *t0,*pt;
face_t f;
cell_t c;
real feta_d,vel_Val[2],cont_Line_Vel,VOF_Normal[2],cap_Num;
real static_Con_Rad,x_Bottom,x_Top,x_Bisect,hoff_Old,hoff_Cur,hoff_New;
real finish_Cond,inv_Hoff=0.0;
int notConverged=MYTRUE;
int itNum=0;
x_Bottom=0.001;
x_Top=2.0;
static_Con_Rad=((static_Con_Ang*M_PI)/180);
while(notConverged)
{
itNum++;
hoff_Old=(Hoff(x_Bottom)-static_Con_Rad);
hoff_Cur=(Hoff(x_Top)-static_Con_Rad);
x_Bisect=(x_Bottom+x_Top)/2.0;
hoff_New=(Hoff(x_Bisect)-static_Con_Rad);
finish_Cond=fabs(1-(x_Bisect/x_Top));
if(finish_Cond<0.00000001 || itNum>10000000)
{
inv_Hoff=x_Bisect;
notConverged=MYFALSE;
}
if((hoff_Old*hoff_New)<0.0)
{
x_Top=x_Bisect;
}
if((hoff_Cur*hoff_New)<=0.0)
{
x_Bottom=x_Bisect;
}
}
begin_f_loop(f,thread)
{
c=F_C0(f,thread);
t0=THREAD_T0(thread);
pt=THREAD_SUB_THREAD(t0,1);
if(C_VOF(c,pt)!=0.0)
{
vel_Val[0]=C_U(c,t0);
vel_Val[1]=C_V(c,t0);
VOF_Normal[0]=C_UDMI(c,t0,0);
VOF_Normal[1]=C_UDMI(c,t0,1);
cont_Line_Vel=NV_DOT(vel_Val,VOF_Normal);
cap_Num=fabs((VISCOSITY*cont_Line_Vel)/SURF_TENS);
if(cap_Num+inv_Hoff<0.0)
{
cap_Num=inv_Hoff;
}
feta_d=((Hoff(cap_Num+inv_Hoff))*180)/M_PI;
F_PROFILE(f,thread,position)=feta_d;
}
else
{
F_PROFILE(f,thread,position)=static_Con_Ang;
}
}
end_f_loop(f,thread)
