24小时热门版块排行榜

返回列表

gky227227

新虫 (正式写手)

应助: 0 (幼儿园)
金币: 10
散金: 81
红花: 3
帖子: 304
在线: 90小时
虫号: 807273
注册: 2009-07-11
专业: 化学反应工程

[求助] k-e方程添加源项报错

大家过年好，向大家求助一个问题，最近想在UDF中给k-e方程添加源项，编译通过了，可是计算时第一步就报错了，FLUENT received fatal signal (ACCESS_VIOLATION)，如图所示，程序的代码如下，还请大家给看一下是什么问题，谢谢。另外，我添加源项的方法是在cell zone condition 中添加的。
DEFINE_SOURCE(Source_kediss,cell,thread,dS,eqn)
{
Thread *mix_thread, *thread_gas, *thread_liq;
cell_t c_back, c_front;
int np,i, j, k, j_temp;
float TEl_g, EDl_g, epsilon, delta_e, lumpted, lumptis, gasfraction, Db_Sauter, Sediss;
float tempb,tempc, timelifej, timelifek, ubtj, ubtk, hbj, hbk, hbjk, lbjk, chl;
float k_bl=500.0, k_b=6.0, alphac=0.8;
float Source_ediss;
FILE *fp;

float FD, vof_gas, y_vel_gas, y_vel_liq, vel_slip, tao, Sk;
float x_vel_liq, x_vel_gas, y_vel_slip, x_vel_slip, ke, ediss;
float sumf=0,temp_x=0, f_lar=0.0, kb_lar=0.0, FD_total=0.0;
float Re[M-1], Eo[M-1], CD[M-1], temp_a[M-1], temp_b[M-1];
float f[M], Num[M];

thread_liq= THREAD_SUB_THREAD(mix_thread, 0);
thread_gas= THREAD_SUB_THREAD(mix_thread, 1);
C_VOF(cell, thread_gas)=MAX(C_VOF(cell, thread_gas), 1.0e-6);
C_VOF(cell, thread_gas)=MIN(C_VOF(cell, thread_gas), 0.64);    /*这个处理如果有自由界面需要去掉*/
vof_gas=C_VOF(cell, thread_gas);
x_vel_liq=C_U(cell, thread_liq);
x_vel_gas=C_U(cell, thread_gas);
y_vel_liq=C_V(cell, thread_liq);
y_vel_gas=C_V(cell, thread_gas);
y_vel_slip=fabs(y_vel_gas-y_vel_liq);
x_vel_slip=fabs(x_vel_gas-x_vel_liq);
vel_slip=pow(y_vel_slip*y_vel_slip+x_vel_slip*x_vel_slip,0.5);

for(i=0;i<M-1;i++)
temp += C_UDSI(cell, thread_gas, i)/db;
Db_Sauter=1/temp;

for(i=0; i<29; i++)
{
C_UDSI(cell, thread_gas, i)=MAX(C_UDSI(cell, thread_gas, i), 0.0);
}


/* set sum(f) to be 1*/
sumf=0.0;
for(i=0; i<29; i++)
{
f=C_UDSI(cell, thread_gas, i);
Num=f*C_VOF(cell, thread_gas)/mb;
sumf += f;
}
for(i=0; i<29; i++)
{
f=f/sumf;
C_UDSI(cell, thread_gas, i)=f;
Num=f*C_VOF(cell, thread_gas)/mb;
}

for(i=0; i<M-1; i++)
{
Eo=9.81*(roul-roug)*pow(db,2.0)/sigma;
      Re=roul*db*vel_slip/visl;
temp_a=24.0/Re*(1.0+0.15*pow(Re,0.687));
temp_b=8.0/3.0*Eo/(Eo+4.0);
      CD=MAX(temp_a,temp_b);
FD_total += C_UDSI(cell, thread_gas, i)*CD/db;
if(db>0.0072)
f_lar += C_UDSI(cell, thread_gas, i);
}

temp_x=50.0*vof_gas*f_lar;
kb_lar=MAX(1.0,temp_x);
FD=vof_gas*roul*3.0/4.0*FD_total*vel_slip/kb_lar;
Sk=FD*vel_slip;

ke= C_K(cell, thread);
ediss= C_D(cell, thread);

tao=pow((Db_Sauter*Db_Sauter/ediss),1/3);
/*tao=Db_Sauter/pow(ke, 0.5);*/

Sediss=2*Sk/tao;
return Sediss;
}

未命名.jpg