关于物性的UDF

tonyhc

我用UDF编实际气体密度方程（BWRS），但在FLUENT不能计算，在UDF中编译好像是通过了，请大家帮我看看，由于我刚学UDF，可能让大家见笑了。
在这里用正割法求密度，直到达到一定的收敛指标，迭代结束。
#include "udf.h"
#define RGASU UNIVERSAL_GAS_CONSTANT
#define A0 118.327853
#define B0 0.04034
#define C0 625766.981722
#define r 0.004347
#define b 0.004391
#define a 4.345709
#define e 0.00005
#define c 74572.151916
#define D0 7024090.36446
#define d 96.134637
#define E0 138123501.518842
#define R 8.3143
#define MWT 17.9811
DEFINE_PROPERTY(cell_density, cell, thread)
{
  int i;
  real temp = C_T(cell, thread);
  real press=C_P(cell, thread);
  real density[10],p4[10],F[10],p1,p2,p3;
  density[1]=0;
  density[2]=press/(R*temp);
  p1=B0*R*temp-A0-C0/(temp*temp)+D0/(temp*temp*temp)-E0/(temp*temp*temp*temp);
  p2=b*R*temp-a-d/temp;
  p3=e*(a+d/temp);
  for(i=1;i<3;i++)
  {p4=c*pow(density,3)*(1+r*density*density)*exp(-r*density*density)/(temp*temp);
  F=density*R*temp+p1*pow(density,2)+p2*pow(density,3)+p3*pow(density,6)+p4;
  }
for(i=3;;i++)
{density=(density[i-2]*F[i-1]-density[i-1]*F[i-2])/(F[i-1]-F[i-2]);
  p4=c*pow(density,3)*(1+r*density*density)*exp(-r*density*density)/(temp*temp);
  F=density*R*temp+p1*pow(density,2)+p2*pow(density,3)+p3*pow(density,6)+p4;
if(fabs(density-density[i-1])<1e-6)
  break;
  }
return (density[i-2]*F[i-1]-density[i-1]*F[i-2])/(F[i-1]-F[i-2]);
}