Error when using k-e model !?

Mika

Hi, I'm simulating with Fluent the development of a flow initially at rest under EM forces.

I wrote an UDF which calculates the distribution of EM forces.

However, I have some difficulties to make the solution converge to a "physical solution": I change the solver parameters, the mesh... but nothing seems to change the numerical solution.

Finally, I decided to run the solution wih a k-e model instead of the laminar one. But, before the end of the first iteration, and only when using the k-e model, I have this error message:

Error: > (greater-than): invalid argument [2]: wrong type [not a number] Error Object: nanq

I verified that my UDF gives always correct numbers, and it is.

Thanks for ideas or solutions...

Mika

PS: here is my UDF, as asked in a previous post.

double zmaxi; double pas;

int NbreAimants; double *TabAimants;

int NbreElectrodes; double *TabElectrodes;

real x_fluent[3]; double X[3]; double F[3];

char path[100]="/legi44p/tourn/perrin/DEA/"; /*char path[100]="/legipc65/home/perrin/DEA/";*/ char path2[100];

FILE *debug_x; /*FILE *debug_y;*/ /* FILE *count;*/

DEFINE_SOURCE(emf_x,cell,thread,dS,eqn) {

int i;

double z;

double force[3];

double j;

int it;

/* permet de recuperer le nombre d'iteration */

it = (nres == 0) ? (0) : ((int) count2[nres - 1]);

if(it<1)

{

strcpy(path2,path);

if((debug_x=fopen(strcat(path2,"fluent/udf/emf7/src/debug_x"),"a"))==NULL)

{

fprintf(stderr,"Erreur a l'ouverture du fichier debug");

}

/* strcpy(path2,path);

if ((count=fopen(strcat(path2,"fluent/udf/emf7/src/count"),"r"))==NULL) {

fprintf(stderr,"Erreur a l'ouverture du fichier count 1");

}

fscanf(count,"%d",&compteur);

fclose(count);

strcpy(path2,path);

if ((count=fopen(strcat(path2,"fluent/udf/emf7/src/count"),"w"))==NULL) {

fprintf(stderr,"Erreur a l'ouverture du fichier count 2");

}

fprintf(count,"%d",++compteur);

fclose(count);*/

lect_data();

lect_aimants();

lect_electrodes();

/* récupération des points du calcul */

C_CENTROID(x_fluent,cell,thread);

X[0]=x_fluent[0]*1000.0;

X[1]=x_fluent[1]*1000.0;

X[2]=x_fluent[2]*1000.0;

z=zmini;

i=0;

force[0]=0;

force[1]=0;

force[2]=0;

while(z<zmaxi)

{

X[2]=z;

calcul_force(X);

force[0]+=F[0];

force[1]+=F[1];

force[2]+=F[2];

z+=pas;

i++;

}

j=(double)i;

F[0]=force[0]/j;

F[1]=force[1]/j;

F[2]=force[2]/j;

fprintf(debug_x,"\nX: %lf, Y: %lf, Z:%lf, force_x: %lf Fx=%f force_y: %lf FY=%fj=%lf\n",X[0],X[1],X[2],force[0],(real)F[0],force[1],(real)F[1],j);

free(TabAimants);

free(TabElectrodes);

/* sauvegarde des valeurs dans une variable */

C_UDMI(cell,thread,0)=(real)F[0];

}

else

{

F[0]=C_UDMI(cell,thread,0);

}

fclose(debug_x);

return (real)F[0]; }

void lect_data() {

FILE *fdata;

char str[25];

strcpy(path2,path);

if ((fdata=fopen(strcat(path2,"fluent/udf/emf7/src/data.dat"),"r")) ==NULL)

{

fprintf(stderr,"Erreur a l'ouverture du fichier data.dat");

}

fscanf(fdata,"%s %lf",str,&conductivite);

pi=M_PI;

fscanf(fdata,"%s %lf",str,&unite);

fscanf(fdata,"%s %lf",str,&zmini);

fscanf(fdata,"%s %lf",str,&zmaxi);

fscanf(fdata,"%s %lf",str,&pas);

fclose(fdata); }

void lect_aimants() {

FILE *faimants;

int i;

int j;

strcpy(path2,path);

if((faimants=fopen(strcat(path2,"fluent/udf/emf7/src/aimants.dat"),"r"))==NULL)

{

fprintf(stderr,"Erreur a l'ouverture du fichier aimants.dat");

}

fscanf(faimants,"%d",&NbreAimants);

TabAimants=malloc(NbreAimants*6*sizeof(double));

for(i=0;i<NbreAimants;i++)

{

for(j=0;j<6;j++)

{

fscanf(faimants,"%lf",(TabAimants+i*6+j));

}

}

fclose(faimants); }

void lect_electrodes() {

FILE *felectrodes;

int i;

int j;

strcpy(path2,path);

if((felectrodes=fopen(strcat(path2,"fluent/udf/emf7/src/electrodes.dat"),"r"))==NULL) {

fprintf(stderr,"Erreur a l'ouverture du fichier electrodes.dat");

}

fscanf(felectrodes,"%d",&NbreElectrodes);

TabElectrodes=malloc(NbreElectrodes*6*sizeof(doubl e));

for(i=0;i<NbreElectrodes;i++)

{

for(j=0;j<6;j++)

{

fscanf(felectrodes,"%lf",(TabElectrodes+i*6+j));

}

}

fclose(felectrodes); } void calcul_force(double *x) {

double B[3];

double E[3];

int i;

calcul_B(x,B);

calcul_E(x,E);

prod_vect(E,B,F);

for(i=0;i<3;i++)

{

F[i]*=conductivite;

} }

void prod_vect(double *A, double *B, double *C) {

*C=(*(A+1)**(B+2)-*(A+2)**(B+1));

*(C+1)=-(*A**(B+2)-*B**(A+2));

*(C+2)=(*A**(B+1)-*B**(A+1)); }

void calcul_B(double *x, double *B) {

int i;

int j;

double k;

double surf[6];

double tmpB[3];

for(j=0;j<3;j++)

{

*(B+j)=0;

*(tmpB+j)=0;

}

for(i=0;i<NbreAimants;i++)

{

for(j=0;j<6;j++)

{

surf[j]=*(TabAimants+i*6+j);

}

k=surf[5];

yonnet(tmpB,x,surf,k);

for(j=0;j<3;j++)

{

B[j]+=tmpB[j];

}

} }

void calcul_E(double *x, double *E) {

int i;

int j;

double k;

double surf[6];

double tmpE[3];

for(j=0;j<3;j++)

{

*(E+j)=0;

*(tmpE+j)=0;

}

for(i=0;i<NbreElectrodes;i++)

{

for(j=0;j<6;j++)

{

surf[j]=*(TabElectrodes+i*6+j);

}

k=surf[5]/(surf[4]*surf[3]*(pow(unite,(double)2))*conductivite);

yonnet(tmpE,x,surf,k);

for(j=0;j<3;j++)

{

E[j]+=tmpE[j];

}

} }

void yonnet(double *C, double *x, double *surf, double k) {

int i;

int j;

double xs[3];

double a;

double b;

double coef;

double S[2];

double T[2];

double R[2][2];

double tmp[3];

for(i=0;i<3;i++)

{

xs[i]=(*(x+i)-*(surf+i))*unite;

}

a=*(surf+3)/2*unite;

b=*(surf+4)/2*unite;

coef=k/(2*pi);

for(i=0;i<2;i++)

{

S[i]=xs[2]-(pow((double)-1,(double)(i+1)))*b;

/* Attention les axes dans le yonnet et dans notre etude ne sont pas les mêmes */

T[i]=xs[0]-(pow((double)-1,(double)(i+1)))*a;

}

for(i=0;i<2;i++)

{

for(j=0;j<2;j++)

{

R[i][j]=(sqrt(pow(S[i],(double)2)+pow(T[j],(double)2)+pow(xs[1],(double)2)));

}

}

for(i=0;i<3;i++)

{

tmp[i]=0;

}

for(i=0;i<2;i++)

{

for(j=0;j<2;j++)

{

tmp[0]+=(pow((double)-1,(double)(i+j))*log((R[i][j]-S[i])));

tmp[1]+=(pow((double)-1,(double)(i+j))*atan((S[i]*T[j]/(R[i][j]*xs[1]))));

tmp[2]+=(pow((double)-1,(double)(i+j))*log((R[i][j]-T[j])));

}

}

for(i=0;i<3;i++)

{

*(C+i)=coef*tmp[i];

}

}