|
|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?注册
x
各位大侠,小弟编写了一个UDF,使用编译型连接。在计算过程中无法识别循环过程,求教如何解决啊?
#include"udf.h"
real G,p1,p2;
DEFINE_ADJUST(inlet,d)
{
real sum_M_A=0.0;
int i,j,k;
real q[200],F[200],F1[200],F2[200],Y[200],Fmax;
real pj[200],dpj[200];
int Ck[3][6]={0,0,0,0,0,0,0,1,0,-1,1,-1,0,0,1,0,-1,1};
real R[6]={0,0.0003824,0.0008637,0.0002366,0.0002366,0.000238475};
real G,Q[6],fjts[6]={0,1,0,0,0,0};
Thread *thread_out;
face_t f;
Domain *domain;
domain=Get_Domain(1);
thread_out=Lookup_Thread(domain,5);
begin_f_loop(f,thread_out)
{
sum_M_A+=F_FLUX(f,thread_out);
}
end_f_loop(f,thread_out)
G=-sum_M_A;
Q[1]=G/1.2;
Q[2]=15;
Q[3]=300+15-G/1.2;
Q[4]=G/1.2+15;
Q[5]=300-G/1.2;
for(i=1;i<=5;i++)
{
pj=fjts*(-1.2*23*0.8*fabs(Q)/80/80+1.2*23*0.8*23/80);
dpj=fjts*(-1.2*23*0.8/80/80);
}
do{
for(i=1;i<=2;i++)
{
F=0;
Y=0;
for(j=1;j<=5;j++)
{
F=F+Ck[j]*R[j]*Q[j]*fabs(Q[j])-Ck[j]*pj[j];
Y=Y+2*Ck[j]*Ck[j]*R[j]*fabs(Q[j])-Ck[j]*Ck[j]*dpj[j];
}
q=-F/Y;
for(j=1;j<=5;j++)
{
if(Ck[j]>0)
Q[j]=Q[j]+q;
else if(Ck[j]<0)
Q[j]=Q[j]-q;
}
F1=0;
F2=0;
for(j=1;j<=5;j++)
{
if(Ck[j]/Q[j]>0)
F1=F1+Ck[j]*R[j]*Q[j]*fabs(Q[j]);
}
for(j=1;j<=5;j++)
{
if(Ck[j]/Q[j]<0)
F2=F2+Ck[j]*R[j]*Q[j]*fabs(Q[j]);
}
for(j=1;j<=5;j++)
{
if(Ck[j]*pj[j]<0)
F1=F1-Ck[j]*pj[j];
}
for(j=1;j<=5;j++)
{
if(Ck[j]*pj[j]>0)
F2=F2-Ck[j]*pj[j];
}
for(j=1;j<=5;j++)
{
pj=fjts*(-1.2*23*0.8*fabs(Q)/80/80+1.2*23*0.8*23/80);
dpj=fjts*(-1.2*23*0.8/80/80);
}
}
Fmax=0;
for(k=1;k<=2;k++)
{
if(fabs(fabs(F1[k]/F2[k])-1)>Fmax)
Fmax=fabs(fabs(F1[k]/F2[k])-1);
}
}
while(Fmax>0.0001);
p1=0.000238475*Q[5]*Q[5]+0.0002366*Q[3]*Q[3]-0.0002366*Q[4]*Q[4]+pj[1]-0.0003824*7/20*Q[1]*Q[1]+0.5*1.2*Q[1]*Q[1]/80/80;
p2=0.0003824*Q[1]*Q[1]/2;
}
DEFINE_PROFILE(in,thread,position)
{
face_t f;
begin_f_loop(f,thread)
{
if(G<=0)
F_PROFILE(f,thread,position)=1;
else
F_PROFILE(f,thread,position)=p1;
}
end_f_loop(f,thread)
}
DEFINE_PROFILE(out,thread,position)
{
face_t f;
begin_f_loop(f,thread)
{
if(G<=0)
F_PROFILE(f,thread,position)=1;
else
F_PROFILE(f,thread,position)=p2;
}
end_f_loop(f,thread)
} |
|