fluent中pbm破碎核udf求解

[CFD_xiaobai]

#include "udf.h"

#include "sg_pb.h"

#include "sg_mphase.h"

#define rhoC 998.2

#define sigma 0.072

#define rho_d 1.2

DEFINE_PB_BREAK_UP_RATE_FREQ(break_up_freq_tav, cell, thread, d_1)
{

  Thread *mixture=THREAD_SUPER_THREAD(thread);

  real k,epsi,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,omega,alpha,psi,itea;

  real pi=3.14159,hf=1.0e-8,h0=1.0e-4;
real C1 = 0.675, C2 = 0.1328;
int i,s=0; 
  Thread *tm;
omega=C_O(cell,mixture);
k=C_K(cell, mixture);
epsi=0.09*k*omega;

alpha=C_VOF(cell,thread);
itea=11.4*pow(10,-4.5)/pow(epsi,1/4);
 
psi=pow(epsi,2/3)*pow(d_1,5/3);

f1=0.1*pow((1+0.2),2)/pow(0.2,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.2,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.2,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.2,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.2,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.2,11/3))                            );
f2=0.1*pow((1+0.3),2)/pow(0.3,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.3,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.3,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.3,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.3,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.3,11/3))                            );
f3=0.1*pow((1+0.4),2)/pow(0.4,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.4,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.4,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.4,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.4,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.4,11/3))                            );
f4=0.1*pow((1+0.5),2)/pow(0.5,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.5,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.5,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.5,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.5,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.5,11/3))                            );
f5=0.1*pow((1+0.6),2)/pow(0.6,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.6,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.6,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.6,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.6,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.6,11/3))                            );
f6=0.1*pow((1+0.7),2)/pow(0.7,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.7,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.7,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.7,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.7,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.7,11/3))                            );
f7=0.1*pow((1+0.8),2)/pow(0.8,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.8,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.8,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.8,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.8,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.8,11/3))                            );
f8=0.1*pow((1+0.9),2)/pow(0.9,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.9,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.9,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.9,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.9,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.9,11/3))                            );
f9=0.2*(1+exp(-0.26*0.00042/psi)+2*exp(-0.1476*0.00042/psi)+2*exp(-0.2038*0.00042/psi)+2*exp(-0.2365*0.00042/psi)+2*exp(-0.2543*0.00042/psi)                            );
f10=(0.05-itea/2)*pow((1+itea),2)/pow(itea,11/3)*(1+exp(-0.26*0.00042/psi/pow(itea,11/3))+2*exp(-0.1476*0.00042/psi/pow(itea,11/3))+2*exp(-0.2038*0.00042/psi/pow(itea,11/3))+2*exp(-0.2365*0.00042/psi/pow(itea,11/3))+2*exp(-0.2543*0.00042/psi/pow(itea,11/3))                            );
f11=(0.1-itea/2)*pow((1+0.1),2)/pow(0.1,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.1,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.1,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.1,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.1,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.1,11/3))                            );





     return 0.923*pow(epsi,1/3)/pow(d_1,2/3)*(1-alpha)*(f1+f2+f3+f4+f5+f6+f7+f8+f9+f10+f11);
}
各位大佬，请问我这个根据罗和安教授的破碎核函数写的破碎频率udf，可以正常导入fluent，计算方法使用fluent自带的CM方法，但是除了bin0，剩下的bin方程残差都一样，求解答。下面是子尺寸分布函数
#include "udf.h"
 #include "sg_pb.h"
 #include "sg_mphase.h"
#define max(a, b) (((a) > (b)) ? (a) : (b))
#define min(a,b) ( ((a)>(b)) ? (b):(a) )
 DEFINE_PB_BREAK_UP_RATE_PDF(break_up_pdf_par, cell, thread, d_1, thread_2, d_2)
 {
Thread *mixture=THREAD_SUPER_THREAD(thread);
     real pdf;
     real f_v = min(d_2,pow(pow(d_1,3)-pow(d_2,3),1/3))/d_1;
     real f_bv =pow(f_v,3.) ;
 real c_f =pow(f_bv,2/3)+pow((1-f_bv),2/3)-1 ;
     real k,epsi,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18,f19,f20,f21,f22,f23,f24,omega,alpha,psi,itea;

  Thread *tm;

omega=C_O(cell,mixture);
k=C_K(cell, mixture);
alpha=C_VOF(cell,thread);
  epsi=0.09*k*omega;
itea=11.4*pow(10,-4.5)/pow(epsi,1/4);
psi=pow(epsi,2/3)*pow(d_1,5/3);



f1=0.1*pow((1+0.2),2)/pow(0.2,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.2,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.2,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.2,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.2,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.2,11/3))                            );
f2=0.1*pow((1+0.3),2)/pow(0.3,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.3,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.3,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.3,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.3,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.3,11/3))                            );
f3=0.1*pow((1+0.4),2)/pow(0.4,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.4,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.4,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.4,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.4,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.4,11/3))                            );
f4=0.1*pow((1+0.5),2)/pow(0.5,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.5,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.5,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.5,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.5,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.5,11/3))                            );
f5=0.1*pow((1+0.6),2)/pow(0.6,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.6,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.6,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.6,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.6,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.6,11/3))                            );
f6=0.1*pow((1+0.7),2)/pow(0.7,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.7,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.7,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.7,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.7,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.7,11/3))                            );
f7=0.1*pow((1+0.8),2)/pow(0.8,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.8,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.8,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.8,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.8,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.8,11/3))                            );
f8=0.1*pow((1+0.9),2)/pow(0.9,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.9,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.9,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.9,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.9,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.9,11/3))                            );
f9=0.2*(1+exp(-0.26*0.00042/psi)+2*exp(-0.1476*0.00042/psi)+2*exp(-0.2038*0.00042/psi)+2*exp(-0.2365*0.00042/psi)+2*exp(-0.2543*0.00042/psi)                            );
f10=(0.05-itea/2)*pow((1+itea),2)/pow(itea,11/3)*(1+exp(-0.26*0.00042/psi/pow(itea,11/3))+2*exp(-0.1476*0.00042/psi/pow(itea,11/3))+2*exp(-0.2038*0.00042/psi/pow(itea,11/3))+2*exp(-0.2365*0.00042/psi/pow(itea,11/3))+2*exp(-0.2543*0.00042/psi/pow(itea,11/3))                            );
f11=(0.1-itea/2)*pow((1+0.1),2)/pow(0.1,11/3)*(1+exp(-0.26*0.00042/psi/pow(0.1,11/3))+2*exp(-0.1476*0.00042/psi/pow(0.1,11/3))+2*exp(-0.2038*0.00042/psi/pow(0.1,11/3))+2*exp(-0.2365*0.00042/psi/pow(0.1,11/3))+2*exp(-0.2543*0.00042/psi/pow(0.1,11/3))                            );
f23=f1+f2+f3+f4+f5+f6+f7+f2+f8+f9+f10+f11;

f12=0.1*pow((1+0.2),2)/pow(0.2,11/3)*exp(-0.00042*c_f/psi/pow(0.2,11/3));
f13=0.1*pow((1+0.3),2)/pow(0.3,11/3)*exp(-0.00042*c_f/psi/pow(0.3,11/3));
f14=0.1*pow((1+0.4),2)/pow(0.4,11/3)*exp(-0.00042*c_f/psi/pow(0.4,11/3));
f15=0.1*pow((1+0.5),2)/pow(0.5,11/3)*exp(-0.00042*c_f/psi/pow(0.5,11/3));
f16=0.1*pow((1+0.6),2)/pow(0.6,11/3)*exp(-0.00042*c_f/psi/pow(0.6,11/3));
f17=0.1*pow((1+0.7),2)/pow(0.7,11/3)*exp(-0.00042*c_f/psi/pow(0.7,11/3));
f18=0.1*pow((1+0.8),2)/pow(0.8,11/3)*exp(-0.00042*c_f/psi/pow(0.8,11/3));
f19=0.1*pow((1+0.9),2)/pow(0.9,11/3)*exp(-0.00042*c_f/psi/pow(0.9,11/3));
f20=0.2*exp(-0.00042*c_f/psi);
f21=(0.05-itea/2)*pow((1+itea),2)/pow(itea,11/3)*exp(-0.00042*c_f/psi/pow(itea,11/3));
f22=(0.1-itea/2)*pow((1+0.1),2)/pow(0.1,11/3)*exp(-0.00042*c_f/psi/pow(0.1,11/3));
f24=f12+f13+f14+f15+f16+f17+f18+f19+f20+f21+f22;



pdf =f24/f23;
     return pdf;
 }

[李东岳]

感谢分享，1024！

[CFD_xiaobai]

@李东岳 请问东岳老师遇到过，class method中除bin0方程外，其余bin方程残差都一致的情况吗？