马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
我用vof法计算气液两相流,想先找出两相界面,并求界面电密度,用DEFINE_ADJUST和DEFINE_ON_DEMAND都试过,能顺利load进去,但一hook这个界面电密度文件fluent就提示致命错误。下面是我的部分程序,请高手指导一下,先谢了。
为方便,附上附件,最后一个udf DEFINE_ADJUST(calm_rhos, domain)出错。
vof_interface_force.txt
(3.23 KB, 下载次数: 0)
28 秒前 上传
点击文件名下载附件
#include "stdio.h"
#include "udf.h"
#include "sg.h"
#include "sg_mphase.h"
#include "math.h"
#include "flow.h"
#include "mem.h"
#include "metric.h"
#include "id.h"
#define EPS_1 4.25 /* 相对介电常数*/
#define EPS_0 8.8541878e-12 /* electrical permittivity 真空介电常数c^2/Nm^2*/
/* Define which user-defined scalars to use. */
enum /*枚举类型名*/
{
phi, /* 电势*/
elrho, /* 电密度*/
N_REQUIRED_UDSI /*存取用户自定义标量输运方程的数目*/
};
enum /*枚举类型名*/
{
Ex,
Ey,
Ez,
vfx,
vfy,
vfz,
rhos, /* 界面电密度*/
N_REQUIRED_UDM /*存取用户自定义标量输运方程的数目*/
};
DEFINE_ADJUST(store_VOF_gradient, d) /*存储vof体积分数的梯度*/
{
Thread *t;
Thread *ppt;
Thread **pt;
cell_t c;
int phase_domain_index=1;
Domain *pDomain = DOMAIN_SUB_DOMAIN(d,phase_domain_index);
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);
mp_thread_loop_c(t,d,pt)
{
if (FLUID_THREAD_P(t))
{
ppt = pt[phase_domain_index];
begin_c_loop (c,t)
{
C_UDMI(c,t,vfx) = C_VOF_G(c,ppt)[0];
C_UDMI(c,t,vfy) = C_VOF_G(c,ppt)[1];
C_UDMI(c,t,vfz) = C_VOF_G(c,ppt)[2];
}
end_c_loop (c,t)
}
}
Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
}
DEFINE_ADJUST(Calc_E_xyz, domain) /*存储电场强度*/
{
Thread *t;
cell_t c;
face_t f;
Thread *t0;
cell_t c0;
real dr0[ND_ND],dr1[ND_ND], xf[ND_ND], dy;
if (n_uds < N_REQUIRED_UDSI)
Internal_Error("not enough user-defined scalars allocated");
if(n_udm < N_REQUIRED_UDM)
Error("Not enough udm allocated\n");
/* Do nothing if gradient isn't allocated yet. */
if (!Data_Valid_P()) return;
thread_loop_c(t,domain)
{
if (FLUID_THREAD_P(t))
if (NULL != THREAD_STORAGE(t,SV_UDS_I(phi)) && NULL != T_STORAGE_R_NV (t,SV_UDSI_G(phi)))
{
begin_c_loop(c,t) /*体循环*/
{
C_UDMI(c,t,Ex) = -1.*C_UDSI_G(c,t,phi)[0]; /*Ex*/
C_UDMI(c,t,Ey) = -1.*C_UDSI_G(c,t,phi)[1];
C_UDMI(c,t,Ez) = -1.*C_UDSI_G(c,t,phi)[2];
}
end_c_loop(c,t)
/* Message("x electric strength: %g\n", C_UDMI(c,t,Ex)); */
}
}
}
DEFINE_ADJUST(calm_rhos, domain) /* 界面电密度,出问题的地方,请帮忙指正*/
{
Thread *t;
Thread *pri_th,*sec_th;
cell_t c;
real xc[ND_ND], EEp[3], EEs[3], nn[3], EE[3], rr;
/* Domain *domain = Get_Domain(1); Get domain pointer */
int ID = 1;
t =Lookup_Thread(domain, ID); /* mixture-level thread pointer */
pri_th = THREAD_SUB_THREAD(t,0);
sec_th = THREAD_SUB_THREAD(t,1);
thread_loop_c(t,domain)
{
begin_c_loop(c,t) /*体循环*/
{
C_CENTROID(xc,c,t);
if (C_VOF(c,sec_th)>0. && C_VOF(c,sec_th)<1.)
{
EEs[0]=C_UDMI(c,sec_th,Ex);/* 次相电场强度 */
EEs[1]=C_UDMI(c,sec_th,Ey);
EEs[2]=C_UDMI(c,sec_th,Ez);
EEp[0]=C_UDMI(c,pri_th,Ex);/* 基相电场强度 */
EEp[1]=C_UDMI(c,pri_th,Ey);
EEp[2]=C_UDMI(c,pri_th,Ez);
nn[0]=C_UDMI(c,sec_th,vfx);/* vof梯度、法向矢量 */
nn[1]=C_UDMI(c,sec_th,vfy);
nn[2]=C_UDMI(c,sec_th,vfz);
NV_V(EEs,=,EEp);
rr=EPS_0*EPS_1;
NV_VS_VS(EE, =,EEp,*,EPS_0, -, EEs,*,rr);
C_UDMI(c,t,rhos) = NV_DOT(nn,EE);
}
else
{
C_UDMI(c,t,rhos) =0;
}
}
end_c_loop(c,t)
}
}
|