REAL GAS UDF
 

I'm modeling a blowdown from a high pressure vessel (7 MPa) to ambient conditions, that go through a nozzle at choked conditions. Does anyone know of a way to model the density accurately, accounting for compressibility? Do you have a UDF? Thanks in advance.

Re: REAL GAS UDF
 

There's a Redlich-Kwong Equation of State based UDF available from Fluent UDF archive, which can be used to model real gases. This may help you...

SRK ?
 

Please wich version of fluent contain the SRK equation of state I am searching it to model h2/o2 combustion ???

Re: SRK ?
 

The UDF is available from the UDF archive at fluentusers.com, and it is a Redlich-Kwong Real Gas UDF, the Soave modified version. However, for H2/O2 combustion I doubt that you will need a real gas model. Have you tried to calculate some compressability factors in order to determine the deviation from an ideal gas at you relevant conditions?

Re: SRK ?
 

Dear Skov,

these are my operating conditions

GH2 P=70 bars , T = 287 K GO2 P=70 bars , T = 100 K

A+

Re: SRK ?
 

What are the respective mole/mass fractions or flow rates -the ratio between H2 and O2.

By the way a real gas model should be available with the coupled solver, thus you do not need to go into UDF modeling of the real gas state...

Re: REAL GAS UDF
 

hello brian, here is the redlich kwong model if you havent found it..

/* The variables below need to be set for a particular gas */

/* CO2 */

/* REAL GAS EQUATION OF STATE MODEL - BASIC VARIABLES */ /* ALL VARIABLES ARE IN SI UNITS! */

#define RGASU UNIVERSAL_GAS_CONSTANT #define PI 3.141592654

#define MWT 44.01 #define PCRIT 7.3834e6 #define TCRIT 304.21 #define ZCRIT 0.2769 #define VCRIT 2.15517e-3 #define NRK 0.77

/* IDEAL GAS SPECIFIC HEAT CURVE FIT */

#define CC1 453.577 #define CC2 1.65014 #define CC3 -1.24814e-3 #define CC4 3.78201e-7 #define CC5 0.00

/* REFERENCE STATE */

#define P_REF 101325 #define T_REF 288.15

Redlich-Kwong Real Gas UDRGM Code Listing

/************************************************** ************/ /* */ /* User-Defined Function: Redlich-Kwong Equation of State */ /* for Real Gas Modeling */ /* */ /* Author: Frank Kelecy */ /* Date: May 2003 */ /* Version: 1.02 */ /* */ /* This implementation is completely general. */ /* Parameters set for CO2. */ /* */ /************************************************** ************/

#include "udf.h" #include "stdio.h" #include "ctype.h" #include "stdarg.h"

/* The variables below need to be set for a particular gas */

/* CO2 */

/* REAL GAS EQUATION OF STATE MODEL - BASIC VARIABLES */ /* ALL VARIABLES ARE IN SI UNITS! */

#define RGASU UNIVERSAL_GAS_CONSTANT #define PI 3.141592654

#define MWT 44.01 #define PCRIT 7.3834e6 #define TCRIT 304.21 #define ZCRIT 0.2769 #define VCRIT 2.15517e-3 #define NRK 0.77

/* IDEAL GAS SPECIFIC HEAT CURVE FIT */

#define CC1 453.577 #define CC2 1.65014 #define CC3 -1.24814e-3 #define CC4 3.78201e-7 #define CC5 0.00

/* REFERENCE STATE */

#define P_REF 101325 #define T_REF 288.15

/* OPTIONAL REFERENCE (OFFSET) VALUES FOR ENTHALPY AND ENTROPY */

#define H_REF 0.0 #define S_REF 0.0

static int (*usersMessage)(char *,...); static void (*usersError)(char *,...);

/* Static variables associated with Redlich-Kwong Model */

static double rgas, a0, b0, c0, bb, cp_int_ref;

DEFINE_ON_DEMAND(I_do_nothing) { /* this is a dummy function to allow us */ /* to use the compiled UDFs utility */ }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_error */ /*------------------------------------------------------------*/

void RKEOS_error(int err, char *f, char *msg) { if (err)

usersError("RKEOS_error (%d) from function: %s\n%s\n",err,f,msg); }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_Setup */ /*------------------------------------------------------------*/

void RKEOS_Setup(Domain *domain, char *filename,

void (*messagefunc)(char *format, ...),

void (*errorfunc)(char *format, ...)) {

rgas = RGASU/MWT;

a0 = 0.42747*rgas*rgas*TCRIT*TCRIT/PCRIT;

b0 = 0.08664*rgas*TCRIT/PCRIT;

c0 = rgas*TCRIT/(PCRIT+a0/(VCRIT*(VCRIT+b0)))+b0-VCRIT;

bb = b0-c0;

cp_int_ref = CC1*log(T_REF)+T_REF*(CC2+

T_REF*(0.5*CC3+T_REF*(0.333333*CC4+0.25*CC5*T_REF) ));

usersMessage = messagefunc;

usersError = errorfunc;

usersMessage("\nLoading Redlich-Kwong Library: %s\n", filename); }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_pressure */ /* Returns density given T and density */ /*------------------------------------------------------------*/

double RKEOS_pressure(double temp, double density) {

double v = 1./density;

double afun = a0*pow(TCRIT/temp,NRK);

return rgas*temp/(v-bb)-afun/(v*(v+b0)); }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_spvol */ /* Returns specific volume given T and P */ /*------------------------------------------------------------*/

double RKEOS_spvol(double temp, double press) {

double a1,a2,a3;

double vv,vv1,vv2,vv3;

double qq,qq3,sqq,rr,tt,dd;

double afun = a0*pow(TCRIT/temp,NRK);

a1 = c0-rgas*temp/press;

a2 = -(bb*b0+rgas*temp*b0/press-afun/press);

a3 = -afun*bb/press;

/* Solve cubic equation for specific volume */

qq = (a1*a1-3.*a2)/9.;

rr = (2*a1*a1*a1-9.*a1*a2+27.*a3)/54.;

qq3 = qq*qq*qq;

dd = qq3-rr*rr;

/* If dd < 0, then we have one real root */

/* If dd >= 0, then we have three roots -> choose largest root */

if (dd < 0.) {

tt = sqrt(-dd)+pow(fabs(rr),0.333333);

vv = (tt+qq/tt)-a1/3.;

} else {

tt = acos(rr/sqrt(qq3));

sqq = sqrt(qq);

vv1 = -2.*sqq*cos(tt/3.)-a1/3.;

vv2 = -2.*sqq*cos((tt+2.*PI)/3.)-a1/3.;

vv3 = -2.*sqq*cos((tt+4.*PI)/3.)-a1/3.;

vv = (vv1 > vv2) ? vv1 : vv2;

vv = (vv > vv3) ? vv : vv3;

}

return vv; }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_density */ /* Returns density given T and P */ /*------------------------------------------------------------*/

double RKEOS_density(double temp, double press, double yi[]) {

return 1./RKEOS_spvol(temp, press); /* (Kg/m3) */ }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_dvdp */ /* Returns dv/dp given T and rho */ /*------------------------------------------------------------*/

double RKEOS_dvdp(double temp, double density) {

double a1,a2,a1p,a2p,a3p,v,press;

double afun = a0*pow(TCRIT/temp,NRK);

press = RKEOS_pressure(temp, density);

v = 1./density;

a1 = c0-rgas*temp/press;

a2 = -(bb*b0+rgas*temp*b0/press-afun/press);

a1p = rgas*temp/(press*press);

a2p = a1p*b0-afun/(press*press);

a3p = afun*bb/(press*press);

return -(a3p+v*(a2p+v*a1p))/(a2+v*(2.*a1+3.*v)); }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_dvdt */ /* Returns dv/dT given T and rho */ /*------------------------------------------------------------*/

double RKEOS_dvdt(double temp, double density) {

double a1,a2,dadt,a1t,a2t,a3t,v,press;

double afun = a0*pow(TCRIT/temp,NRK);

press = RKEOS_pressure(temp, density);

v = 1./density;

dadt = -NRK*afun/temp;

a1 = c0-rgas*temp/press;

a2 = -(bb*b0+rgas*temp*b0/press-afun/press);

a1t = -rgas/press;

a2t = a1t*b0+dadt/press;

a3t = -dadt*bb/press;

return -(a3t+v*(a2t+v*a1t))/(a2+v*(2.*a1+3.*v)); }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_Cp_ideal_gas */ /* Returns ideal gas specific heat given T */ /*------------------------------------------------------------*/

double RKEOS_Cp_ideal_gas(double temp) {

return (CC1+temp*(CC2+temp*(CC3+temp*(CC4+temp*CC5)))); }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_H_ideal_gas */ /* Returns ideal gas specific enthalpy given T */ /*------------------------------------------------------------*/

double RKEOS_H_ideal_gas(double temp) {

return temp*(CC1+temp*(0.5*CC2+temp*(0.333333*CC3+

temp*(0.25*CC4+temp*0.2*CC5)))); }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_specific_heat */ /* Returns specific heat given T and rho */ /*------------------------------------------------------------*/

double RKEOS_specific_heat(double temp, double density, double yi[]) {

double delta_Cp,press,v,dvdt,dadt;

double afun = a0*pow(TCRIT/temp,NRK);

press = RKEOS_pressure(temp, density);

v = 1./density;

dvdt = RKEOS_dvdt(temp, density);

dadt = -NRK*afun/temp;

delta_Cp = press*dvdt-rgas-dadt*(1.+NRK)/b0*log((v+b0)/v)

+ afun*(1.+NRK)*dvdt/(v*(v+b0));

return RKEOS_Cp_ideal_gas(temp)+delta_Cp; /* (J/Kg-K) */ }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_enthalpy */ /* Returns specific enthalpy given T and rho */ /*------------------------------------------------------------*/

double RKEOS_enthalpy(double temp, double density, double yi[]) {

double delta_h ,press, v;

double afun = a0*pow(TCRIT/temp,NRK);

press = RKEOS_pressure(temp, density);

v = 1./density;

delta_h = press*v-rgas*temp-afun*(1+NRK)/b0*log((v+b0)/v);

return H_REF+RKEOS_H_ideal_gas(temp)+delta_h; /* (J/Kg) */ }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_entropy */ /* Returns entropy given T and rho */ /*------------------------------------------------------------*/

double RKEOS_entropy(double temp, double density, double yi[]) {

double delta_s,v,v0,dadt,cp_integral;

double afun = a0*pow(TCRIT/temp,NRK);

cp_integral = CC1*log(temp)+temp*(CC2+temp*(0.5*CC3+

temp*(0.333333*CC4+0.25*CC5*temp)))

- cp_int_ref;

v = 1./density;

v0 = rgas*temp/P_REF;

dadt = -NRK*afun/temp;

delta_s = rgas*log((v-bb)/v0)+dadt/b0*log((v+b0)/v);

return S_REF+cp_integral+delta_s; /* (J/Kg-K) */ }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_mw */ /* Returns molecular weight */ /*------------------------------------------------------------*/

double RKEOS_mw(double yi[]) {

return MWT; /* (Kg/Kmol) */ }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_speed_of_sound */ /* Returns s.o.s given T and rho */ /*------------------------------------------------------------*/

double RKEOS_speed_of_sound(double temp, double density, double yi[]) {

double cp = RKEOS_specific_heat(temp, density, yi);

double dvdt = RKEOS_dvdt(temp, density);

double dvdp = RKEOS_dvdp(temp, density);

double v = 1./density;

double delta_c = -temp*dvdt*dvdt/dvdp;

return sqrt(cp/((delta_c-cp)*dvdp))*v; /* m/s */ }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_rho_t */ /*------------------------------------------------------------*/

double RKEOS_rho_t(double temp, double density, double yi[]) {

return -density*density*RKEOS_dvdt(temp, density); }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_rho_p */ /*------------------------------------------------------------*/

double RKEOS_rho_p(double temp, double density, double yi[]) {

return -density*density*RKEOS_dvdp(temp, density); }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_enthalpy_t */ /*------------------------------------------------------------*/

double RKEOS_enthalpy_t(double temp, double density, double yi[]) {

return RKEOS_specific_heat(temp, density, yi); }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_enthalpy_p */ /*------------------------------------------------------------*/

double RKEOS_enthalpy_p(double temp, double density, double yi[]) {

double v = 1./density;

double dvdt = RKEOS_dvdt(temp, density);

return v-temp*dvdt; }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_viscosity */ /*------------------------------------------------------------*/

double RKEOS_viscosity(double temp, double density, double yi[]) {

double mu,tr,tc,pcatm;

tr = temp/TCRIT;

tc = TCRIT;

pcatm = PCRIT/101325.;

mu = 6.3e-7*sqrt(MWT)*pow(pcatm,0.6666)/pow(tc,0.16666)*

(pow(tr,1.5)/(tr+0.8));

return mu; }

/*------------------------------------------------------------*/ /* FUNCTION: RKEOS_thermal_conductivity */ /*------------------------------------------------------------*/

double RKEOS_thermal_conductivity(double temp, double density, double yi[]) {

double cp, mu;

cp = RKEOS_Cp_ideal_gas(temp);

mu = RKEOS_viscosity(temp, density, yi);

return (cp+1.25*rgas)*mu; }

/* Export Real Gas Functions to Solver */

UDF_EXPORT RGAS_Functions RealGasFunctionList = {

RKEOS_Setup, /* initialize */

RKEOS_density, /* density */

RKEOS_enthalpy, /* enthalpy */

RKEOS_entropy, /* entropy */

RKEOS_specific_heat, /* specific_heat */

RKEOS_mw, /* molecular_weight */

RKEOS_speed_of_sound, /* speed_of_sound */

RKEOS_viscosity, /* viscosity */

RKEOS_thermal_conductivity, /* thermal_conductivity */

RKEOS_rho_t, /* drho/dT |const p */

RKEOS_rho_p, /* drho/dp |const T */

RKEOS_enthalpy_t, /* dh/dT |const p */

RKEOS_enthalpy_p /* dh/dp |const T */ };