这是一个自娱自乐的FLUENT udf分享
自制udf分享(以及寻求帮助)
Warning!!!!!!!!!!!!!!!!!!!!!屎山代码警告,意义不明的变量定义警告
Udf功能:粘性函数1,2通过控制流体粘性略微增加了欧拉多相流计算流体注入容器时的稳定性(出现发散迹象时能救回来,但过于离谱的边界条件该发散还是会发散)
速度边界函数:根据第二相流体体积变化率实时调整速度边界入口速度。
迭代结束执行函数:当某一参数满足设定标准时直接中断计算,这里设置的是当出口质量流大于0.05kg/s时停止计算。
#include "udf.h"
#include "math.h"
#include "sg_udms.h"
#include "unsteady.h"
#include "sg_mphase.h"
#define standard 20
#define inlet_zone 9049
#define mass_rate 0.5
#define acc 0.001
#define dec 0.0001
#define outlet_zone1 19
#define outlet_zone2 20
#define flux_standard 0.05
DEFINE_PROPERTY(pp_viscosity,c,t)
{
real SR=C_STRAIN_RATE_MAG(c,t);
real n=0.511;
real k=1800;
real YSR;
real PV=500;
real YP=300;
real EV;
real MAG=0;
YSR=YP/(k-PV);
EV=1;
MAG=sqrt(pow(C_U(c,t),2)+pow(C_V(c,t),2)+pow(C_W(c,t),2));
YSR=YP/(k-PV);
if (SR<YSR)
{EV=YP*((2*YSR-SR)/(YSR*YSR))+k*((2-n)+(n-1)*(SR/YSR));}
else
{EV=(YP/SR)+k*(exp((n-1)*log(SR/YSR)));}
if(MAG>standard)
{EV=MAG/standard*mu*EV;}
return EV;
}
Pp材料的粘性系数函数,本体为非牛顿流体材料的粘性系数,最后增加了一个判定标准,所有速度大于standard变量的pp材料粘性都会随速度线性增大,standard变量已在函数顶部列出可以手动修改。
DEFINE_PROPERTY(air_viscosity,c,t)
{real EV;
real MAG=0;
EV=0.0017894;
MAG=sqrt(pow(C_U(c,t),2)+pow(C_V(c,t),2)+pow(C_W(c,t),2));
if(MAG>standard)
{EV=MAG/standard*mu*EV*3;}
return EV;
}
Air材料的粘性系数函数,同pp材料,对standard变量的修改同样会对本函数生效。
DEFINE_PROFILE(inlet_v_eular,t,i)
{Domain *d;
Thread *tt;
cell_t c;
face_t f;
real a=0;
real b=0.72;
real v=0.00373;
real h=0;
d=Get_Domain(3);
tt=Lookup_Thread(d,inlet_zone);
begin_c_loop(c,tt)
a+=C_VOLUME(c,tt)*C_VOF(c,tt);
end_c_loop(c,tt)
#if RP_NODE
a=PRF_GRSUM1(a);
#endif
begin_c_loop(c,tt)
if(C_UDMI(c,tt,0)==0)
{C_UDMI(c,tt,0)=a;}
if(C_UDMI(c,tt,1)==0)
{C_UDMI(c,tt,1)=a;}
if(C_UDMI(c,tt,2)<0.72)
{C_UDMI(c,tt,2)=b;}
if(N_TIME>3)
{if(N_TIME%2)
{v=2000*(a-C_UDMI(c,tt,0))/PREVIOUS_TIMESTEP;
v=sqrt(pow(v,2));
if(v<mass_rate&v!=0)
{C_UDMI(c,tt,2)=acc+C_UDMI(c,tt,2);}
if(v>mass_rate)
{C_UDMI(c,tt,2)=-dec+C_UDMI(c,tt,2);}
C_UDMI(c,tt,1)=a;}
if(N_TIME%2==0)
{v=2000*(a-C_UDMI(c,tt,1))/PREVIOUS_TIMESTEP;
v=sqrt(pow(v,2));
if(v<mass_rate&v!=0)
{C_UDMI(c,tt,2)=acc+C_UDMI(c,tt,2);}
if(v> mass_rate)
{C_UDMI(c,tt,2)=-dec+C_UDMI(c,tt,2);}
C_UDMI(c,tt,0)=a;}}
b=C_UDMI(c,tt,2);
end_c_loop(c,tt)
b=PRF_GRHIGH1(b);
begin_f_loop(f,t)
F_PROFILE(f,t,i)=b;
end_f_loop(f,t)}
Inlet_v_eluar函数主要通过监控流场域中pp材料的质量实时变化率修改进出口流速,使质量变化率保持恒定,重要的变量为inlet_zone,mass_rate,acc,dec,这些变量均在顶部中列出,其中在更换不同模型的算例时inlet_zone必须改为对应模型的流场域ID,错误的ID将导致算例崩溃,算例ID可在cell_room_condition选项下,对应域名字后面查看。mass_rate变量为目标质量变化率,如果流场中材料变化率低于mass_rate,则入口速度每次迭代都会加上acc的值,直至变化率等于mass_rate,相反若高于mass_rate,则每次迭代会减去dec的值。
DEFINE_EXECUTE_AT_END(execute_vof_calculate)
{Domain *d;
Thread *t1;
Thread *t2;
cell_t c;
real a=0;
real flux=0;
real flag=0;
face_t f;
real area=0;
real farea[ND_ND];
int outlet_zone1=19;
int zone_ID2=20;
d=Get_Domain(3);
t1=Lookup_Thread(d,outlet_zone1);
t2=Lookup_Thread(d,zone_ID2);
begin_f_loop(f,t1)
F_AREA(farea,f,t1);
area+=NV_MAG(farea);
a+=F_FLUX(f,t1);
end_f_loop(f,t1)
begin_f_loop(f,t2)
F_AREA(farea,f,t2);
area+=NV_MAG(farea);
a+=F_FLUX(f,t2);
end_f_loop(f,t2)
#if RP_NODE
flux=PRF_GRSUM1(a)*1e9/2000*CURRENT_TIMESTEP;
area=PRF_GRSUM1(area)*1e6;
#endif
node_to_host_real_1(flux);
#if RP_HOST
if(flux>flux_standard)
{RP_Set_Integer("interruptflag",9);}
#endif
}
execute_vof_calculate函数为vof法的定量停止函数,需要配合Fluent的command指令和TUI命令一同使用,具体指令会在底部列出。execute_vof_calculate函数的重要变量是outlet_zone1,outlet_zone2和flux_standard,均在顶部列出。该函数具体作用是当单位时间内流出指定边界的体积流量高于flux_standard时,会更改interruptflag变量值将其置为9,然后Fluent内部指令检测到interruptflag变量变化会立即终止计算,从而达到定量停止的目的。当更换不同模型时必须将outlet_zone1和outlet_zone2改为对应出口边界的ID,具体ID可在Fluent中boundary condition选下下出口边界名后部查看。
Command指令如下所示:
(if(>(%rpgetvar 'interruptflag) 0)(set! mstop? #t))
TUI命令如下所示:
(rp-var-define 'interruptflag 0 'interger #f)
其中TUI命令每次启动都需要在console指令框中重新输入。
关于速度入口函数,因为找不到体积关于时间的导数宏从而足足定义了三个udmi位辅助计算,从而使得计算效率奇低,如果哪个大佬有办法改进的话麻烦私我一下,非常感谢,另外如果看的人多的话,或许可能也许。。会出一个实际操作教学?如果懒病不犯的话(笑)