[ANSYS, Inc. Logo] return to home search
next up previous contents index

8.2.6 User-Defined Real Gas Models

This section contains examples of UDFs that can be used to customize user-defined real gas models. See this section in the separate User's Guide for the UDRGM overview, limitations, and details on how to set up, build and load a library of user-defined real gas functions.



UDRGM Example: Redlich-Kwong Equation of State

This section describes another example of a user-defined real gas model. You can use this example as the basis for your own UDRGM code. In this example, the Redlich-Kwong equation of state is used in the UDRGM.

This section summarizes the equations used in developing the UDRGM for the Redlich-Kwong equation of state. The model is based on a modified form of the Redlich-Kwong equation of state described in [ 1]. The equations used in this UDRGM will be listed in the sections below.

The following nomenclature applies to this section:


  $a(T)$ = Redlich-Kwong temperature function
  $c$ = speed of sound
  $C_p$ = specific heat
  $H$ = enthalpy
  $n$ = exponent in function $a(T)$
  $p$ = pressure
  $R$ = universal gas constant/molecular weight
  $T$ = temperature
  $S$ = entropy
  $V$ = specific volume
  $\rho$ = density

The superscript 0 designates a reference state, and the subscript $c$ designates a critical point property.

Specific Volume and Density

The Redlich-Kwong equation of state can be written in the following form:


p = \frac{RT}{(V-\tilde{b})} - \frac{a(T)}{V(V+b_0)} (8.2-12)

where

V = \frac{1}{\rho}, \; \; a(T) = a_0\left(\frac{T_c}{T}\right)^n, \; \; a_0 = 0.42747 \frac{R^2T_{c}^{2}}{p_c},


b_0 = 0.08664 \frac{RT_c}{p_c}, \; \; c_0 = \frac{RT_c}{p_c+\frac{a_0}{V_c(V_c+b_0)}}+b_0-V_c, \; \; \tilde{b} = b_0-c_0

Since the real gas model in ANSYS FLUENT requires a function for density as a function of pressure and temperature, Equation  8.2-12 must be solved for the specific volume (from which the density can be easily obtained). For convenience, Equation  8.2-12 can be written as a cubic equation for specific volume as follows:


V^{3} +a_1V^{2} +a_2V + a_3 = 0 (8.2-13)

where

a_1 = c_0 - \frac{RT}{p}, \;\; a_2 = -\left(\tilde{b}b_0 + \f... ... \frac{a(T)}{p}\right), \;\; a_3 = -\frac{a(T) \; \tilde{b}}{p}

Equation  8.2-13 is solved using a standard algorithm for cubic equations (see [ 10] for details). In the UDRGM code, the cubic solution is coded to minimize the number of floating point operations. This is critical for optimal performance, since this function gets called numerous times during an iteration cycle.

It should be noted that the value of the exponent, $n$, in the function $a(T)$ will depend on the substance. A table of values can be found in [ 1] for some common substances. Alternatively, [ 1] states that values of $n$ are well correlated by the empirical equation


n = 0.4986 + 1.1735\omega + 0.475\omega^{2} (8.2-14)

where $\omega$ is the acentric factor, defined as


\omega = -\log \left(\frac{p_v(T)}{p_c}\right) -1.0 (8.2-15)

In the above equation, $p_v(T)$ is the saturated vapor pressure evaluated at temperature $T=0.7T_c$.

Derivatives of Specific Volume and Density

The derivatives of specific volume with respect to temperature and pressure can be easily determined from Equation  8.2-12 using implicit differentiation. The results are presented below:


$\displaystyle \left(\frac{\partial V}{\partial p}\right)_T$ $\textstyle =$ $\displaystyle -\frac{(a_1)'_p V^{2} + (a_2)'_p V + (a_3)'_p}{3V^{2} + 2a_1V +a_2}$ (8.2-16)
       
$\displaystyle \left(\frac{\partial V}{\partial T}\right)_p$ $\textstyle =$ $\displaystyle -\frac{(a_1)'_T V^{2} + (a_2)'_T V + (a_3)'_T}{3V^{2} + 2a_1V +a_2}$ (8.2-17)

where

(a_1)'_p = \frac{RT}{p^{2}}, \;\; (a_2)'_p = \frac{RTb_0-a(T)}{p^{2}}, \;\; (a_3)'_p = \frac{a(T) \; \tilde{b}}{p^{2}}


(a_1)'_T = -\frac{R}{p}, \;\; (a_2)'_T = \frac{-Rb_0 + \frac{... ... \frac{\tilde{b}}{p}, \;\; \frac{da(T)}{dT} = -n \frac{a(T)}{T}

The derivatives of density can be obtained from the above using the relations


\left(\frac{\partial \rho}{\partial p}\right)_T = - \rho^2\left(\frac{\partial V}{\partial p}\right)_T (8.2-18)


\left(\frac{\partial \rho}{\partial T}\right)_p = - \rho^2\left(\frac{\partial V}{\partial T}\right)_p (8.2-19)

Specific Heat and Enthalpy

Following [ 1], enthalpy for a real gas can be written


H = H^{0}(T) + pV - RT - \frac{a(T)}{b_0}(1+n)\ln \left(\frac{V+b_0}{V}\right) (8.2-20)

where $H^{0}(T)$ is the enthalpy function for a thermally perfect gas (i.e., enthalpy is a function of temperature alone). In the present case, we employ a fourth-order polynomial for the specific heat for a thermally perfect gas [ 8]


C_{p}^{0}(T) = C_1 + C_2T + C_3T^2 + C_4T^3 + C_5T^4 (8.2-21)

and obtain the enthalpy from the basic relation


H^{0}(T) = \int_{T_0}^TC_{p}^0(T)dT (8.2-22)

The result is


H^{0}(T) = C_1T + \frac{1}{2}C_2T^2 + \frac{1}{3}C_3T^3 + \frac{1}{4}C_4T^4 + \frac{1}{5}C_5T^5 - H^{0}(T^{0}) (8.2-23)

Note that $H^{0}(T^{0})$ is the enthalpy at a reference state $(p^{0},T^{0})$, which can be chosen arbitrarily.

The specific heat for the real gas can be obtained by differentiating Equation  8.2-20 with respect to temperature (at constant pressure):


C_p = \left(\frac{\partial H}{\partial T}\right)_p (8.2-24)

The result is


C_p = C_{p}^{0}(T) + p\left(\frac{\partial V}{\partial T}\ri... ...)\frac{\left(\frac{\partial V}{\partial T}\right)_p}{V(V+b_0)} (8.2-25)

Finally, the derivative of enthalpy with respect to pressure (at constant temperature) can be obtained using the following thermodynamic relation [ 8]:


\left(\frac{\partial H}{\partial p}\right)_T = V-T\left(\frac{\partial V}{\partial T}\right)_p (8.2-26)

Entropy

Following [ 1], the entropy can be expressed in the form


S = S^{0}(T,p^{0}) + R \ln \left(\frac{V-\tilde{b}}{V^{0}}\r... ...(\frac{da(T)}{dT}\right)}{b_0}\ln \left(\frac{V+b_0}{V}\right) (8.2-27)

where the superscript $0$ again refers to a reference state where the ideal gas law is applicable. For an ideal gas at a fixed reference pressure, $p^0$ , the entropy is given by


S^{0}(T,p^{0}) = S(T^{0},p^{0}) + \int_{T^{0}}^{T}\frac{C_{p}^{0}(T)}{T}dT (8.2-28)

Note that the pressure term is zero since the entropy is evaluated at the reference pressure. Using the polynomial expression for specific heat, Equation  8.2-21, Equation  8.2-28 becomes


S^{0}(T,p^{0}) = S(T^{0},p^{0}) + C_1\ln(T) + C_2T + \frac{1}{2}C_3T^2 + \frac{1}{3}C_4T^3 + \frac{1}{4}C_5T^4 - f(T^{0}) (8.2-29)

where $f(T^{0})$ is a constant, which can be absorbed into the reference entropy $S(T^0,p^0)$.

Speed of Sound

The speed of sound for a real gas can be determined from the thermodynamic relation


c^2 = \left(\frac{\partial p}{\partial \rho}\right)_S = - \l... ...ght)\frac{V^{2}}{\left(\frac{\partial V}{\partial p}\right)_T} (8.2-30)

Noting that,


C_p - C_V = - T \left( \frac{\partial V}{\partial T} \right\... ...t)^{^{^{^{2}}}} \frac{\partial p}{\partial v }\right\vert _{T} (8.2-31)

we can write the speed of sound as


c = V \sqrt{-\left(\frac{C_p}{C_p-\Delta C}\right)\frac{1}{\left(\frac{\partial V}{\partial p}\right)_T}} (8.2-32)

Viscosity and Thermal Conductivity

The dynamic viscosity of a gas or vapor can be estimated using the following formula from [ 2]:


\mu(T) = 6.3\times10^{7}\frac{M_w^{0.5} p^{0.6666}_c}{T^{0.1666}_c} \left(\frac{T^{1.5}_r}{T_r+0.8}\right) (8.2-33)

Here, $T_r$ is the reduced temperature


T_r = \frac{T}{T_c} (8.2-34)

and $M_w$ is the molecular weight of the gas. This formula neglects the effect of pressure on viscosity, which usually becomes significant only at very high pressures.

Knowing the viscosity, the thermal conductivity can be estimated using the Eucken formula [ 4]:


k = \mu\left(C_p+\frac{5}{4}R\right) (8.2-35)

It should be noted that both Equation  8.2-33 and 8.2-35 are simple relations, and therefore may not provide satisfactory values of viscosity and thermal conductivity for certain applications. You are encouraged to modify these functions in the UDRGM source code if alternate formulae are available for a given gas.

Using the Redlich-Kwong Real Gas UDRGM

Using the Redlich-Kwong Real Gas UDRGM simply requires the modification of the top block of #define macros to provide the appropriate parameters for a given substance. An example listing for CO $_2$ is given below. The parameters required are:


  MWT = Molecular weight of the substance
  PCRIT = Critical pressure (Pa)
  TCRIT = Critical temperature (K)
  ZCRIT = Critical compressibility factor
  VCRIT = Critical specific volume (m $^3$/kg)
  NRK = Exponent of $a(T)$ function
  CC1, CC2, CC3, CC4, CC5 = Coefficients of $C_p(T)$ polynomial curve fit
  P_REF = Reference pressure (Pa)
  T_REF = Reference temperature (K)

The coefficients for the ideal gas specific heat polynomial were obtained from [ 8] (coefficients for other substances are also provided in [ 8]). After the source listing is modified, the UDRGM C code can be recompiled and loaded into ANSYS FLUENT in the manner described earlier. 

/* 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, cxboolean vapor_phase, char *filename, int
                 (*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 P, 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 P, 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 P, 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 P, 
							double yi[])
{
   double cp = RKEOS_specific_heat(temp, density, P, 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 P, double yi[])
{
   return -density*density*RKEOS_dvdt(temp, density);
}

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

double RKEOS_rho_p(double temp, double density, double P, double yi[])
{
   return -density*density*RKEOS_dvdp(temp, density);
}

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

double RKEOS_enthalpy_t(double temp, double density, double P, double yi[])
{
   return RKEOS_specific_heat(temp, density, P, yi);
}


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

double RKEOS_enthalpy_p(double temp, double density, double P, 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 P, 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 P, 
								  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       */
};



UDRGM Example: Multiple-Species Real Gas Model

This is a simple example for multiple-species real gas model that provide you with a template which you can use to write a more complex multiple-species UDRGM.

In this example a fluid material is defined in the setup function as a mixture of four species (H20, N2, O2, CO2). The equation of state was the simple ideal gas equation of state. The other thermodynamic properties where defined by an ideal-gas mixing law.

Other auxiliary functions are written to provide individual species property to the principle function set.

The example also provide numerical method of computing $\frac{d\rho}{dT}$, $\frac{d\rho}{dp}$, $\frac{dh}{dT}$, and $\frac{dh}{dp}$ . 

/* 
 *sccs id: @(#)real_ideal.c	1.10 Copyright 1900/11/09 ANSYS, Inc.
 */

/*
 *      Copyright 1988-1998 ANSYS, Inc.
 *      All Rights Reserved
 *
 *      This is unpublished proprietary source code of ANSYS, Inc.
 *      It is protected by U.S. copyright law as an unpublished work
 *      and is furnished pursuant to a written license agreement.  It
 *      is considered by ANSYS, Inc. to be confidential and may not be
 *      used, copied, or disclosed to others except in accordance with
 *      the terms and conditions of the license agreement.
 */


/*
 *  Windows Warning!!!  Including udf.h is for getting definitions for
 *                  ANSYS FLUENT constructs such as Domain.   You must 
 *                  NOT reference any ANSYS FLUENT globals directly from 
 *                  within this module nor link this against any ANSYS 
 *                  FLUENT libs, doing so will cause dependencies on a  
 *                  specific ANSYS FLUENT binary such as fl551.exe and 
 *                  thus won't be version independent.
 */

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

#if RP_DOUBLE
#define SMALL 1.e-20
#else
#define SMALL 1.e-10
#endif

#define NCMAX 20
#define NSPECIE_NAME 80 


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


static double ref_p, ref_T;
static char   gas[NCMAX][NSPECIE_NAME];
static int    n_specs;

double Mixture_Rgas(double yi[]);
double Mixture_pressure(double Temp, double Rho, double yi[]);
double Mw_i(int i);
double Cp_i(double T, double r, int i);
double K_i(double T, double r, int i);
double Mu_i(double T, double r, int i);
double Rgas_i(double T, double r, int i);
double Gm_i(double T, double r, int i);


DEFINE_ON_DEMAND(I_do_nothing)
{
  /*
    This is a dummy function
    must be included to allow for the use of the
    ANSYS FLUENT UDF compilation utility
  */
}

void Mixture_error(int err, char *f, char *msg)
{
  if (err) 
    usersError("Mixture_error (%d) from function: %s\n%s\n",err,f,msg);
}

/*******************************************************************/
/* Mixture Functions                                               */ 
/* These are the only functions called from ANSYS FLUENT Code      */ 
/*******************************************************************/
void MIXTURE_Setup(Domain *domain, cxboolean vapor_phase, char *specielist,
                   int  (*messagefunc)(char *format, ...),
                   void (*errorfunc)(char *format, ...))
{
  /* This function will be called from ANSYS FLUENT after the
     UDF library has been loaded.
     User must enter the number of species in the mixture
     and the name of the individual species.
  */
  
  int i ;
  usersMessage = messagefunc;
  usersError   = errorfunc;
  ref_p        = ABS_P(RP_Get_Real("reference-pressure"),op_pres);
  ref_T        = RP_Get_Real("reference-temperature");

  if (ref_p == 0.0)
    {
      Message0("\n MIXTURE_Setup: reference-pressure was not set by user \n");
      Message0("\n MIXTURE_Setup: setting reference-pressure to 101325 Pa \n");
      ref_p = 101325.0 ;
    }
  /*====================================================*/
  /*=========   User Input Section =====================*/
  /*====================================================*/
  /*
    Define Number of species & Species name.
    DO NOT use space for naming species
  */
  n_specs = 4 ;     

  (void)strcpy(gas[0],"H2O") ; 
  (void)strcpy(gas[1],"N2")  ; 
  (void)strcpy(gas[2],"O2")  ; 
  (void)strcpy(gas[3],"CO2") ;

  /*====================================================*/
  /*=========   End Of User Input Section ==============*/
  /*====================================================*/

  Message0("\n MIXTURE_Setup: RealGas mixture initialization \n");
  Message0("\n MIXTURE_Setup: Number of Species = %d  \n",n_specs);
  for (i=0; i<n_specs; ++i)
    {      
      Message0("\n MIXTURE_Setup: Specie[%d]         = %s  \n",i,gas[i]);
    }

  /*
    concatenate species name into one string
    and send back to fluent
  */
  strcat(specielist,gas[0]);
  for (i=1; i<n_specs; ++i)
    {
      strcat(specielist," ");
      strcat(specielist,gas[i]);
    }
}

double MIXTURE_density(double Temp, double P, double yi[])
{
  double Rgas =  Mixture_Rgas(yi) ;
  
  double r    =  P/(Rgas*Temp);	/* Density at Temp & P */

  return r;			/* (Kg/m^3) */
}

double MIXTURE_specific_heat(double Temp, double density, double P, 
							 double yi[])
{
  double cp=0.0 ;
  int i ;
    
  for (i=0; i<n_specs; ++i)
    cp += yi[i]*Cp_i(Temp,density,i);
  
  return cp;			/* (J/Kg/K) */
}

double MIXTURE_enthalpy(double Temp, double density, double P, double yi[])
{
  double h=0.0;
  int i ;

  for (i=0; i<n_specs; ++i)
    h  += yi[i]*(Temp*Cp_i(Temp,density,i));
  
  return h;			/* (J/Kg) */
}

double MIXTURE_entropy(double Temp, double density, double P, double yi[])
{
  double s   = 0.0   ;
  double Rgas=0.0 ;

  Rgas = Mixture_Rgas(yi);
  
  s    = MIXTURE_specific_heat(Temp,density,P,yi)*log(Temp/ref_T) -
         Rgas*log(P/ref_p) ;

  return s;			/* (J/Kg/K) */
}

double MIXTURE_mw(double yi[])
{
  double MW, sum=0.0 ;
  int i ;

  for (i=0; i<n_specs; ++i)
    sum += (yi[i]/Mw_i(i)) ;

  MW = 1.0/MAX(sum,SMALL)      ;
  
  return MW;			/* (Kg/Kmol) */
}

double MIXTURE_speed_of_sound(double Temp, double density, double P, 
                              double yi[])
{
  double a, cp, Rgas ;
  
  cp   = MIXTURE_specific_heat(Temp,density,P,yi) ;
  Rgas = Mixture_Rgas(yi) ;


  a    = sqrt(Rgas*Temp* cp/(cp-Rgas) ) ;

  return a ;  /* m/s */
}

double MIXTURE_viscosity(double Temp, double density, double P, double yi[])
{
  double mu=0;
  int i ;

  for (i=0; i<n_specs; ++i)
    mu += yi[i]*Mu_i(Temp,density,i);

  return mu;			/* (Kg/m/s) */
}

double MIXTURE_thermal_conductivity(double Temp, double density, double P, 
                                    double yi[])
{
  double kt=0;
  int i ;

  for (i=0; i<n_specs; ++i)
    kt += yi[i]*K_i(Temp,density,i);
  
  return kt;			/* W/m/K */
}

double MIXTURE_rho_t(double Temp, double density, double P, double yi[])
{

  double drdT   ;	/* derivative of rho w.r.t. Temp */
  double p      ;
  double dT=0.01;

  p     = Mixture_pressure(Temp,density, yi);
  drdT  = (MIXTURE_density(Temp+dT,p,yi) - MIXTURE_density(Temp,p,yi) ) /dT;

  return drdT;			/* (Kg/m^3/K) */
}

double MIXTURE_rho_p(double Temp, double density, double P, double yi[])
{
  double drdp    ;
  double p       ;
  double dp= 5.0 ;
  
  p     = Mixture_pressure(Temp,density, yi);
  drdp = (MIXTURE_density(Temp,p+dp,yi) - MIXTURE_density(Temp,p,yi) ) /dp;
  
  return drdp;			/* (Kg/m^3/Pa) */
}

double MIXTURE_enthalpy_t(double Temp, double density, double P, double yi[])
{
  double dhdT     ;
  double p        ;
  double rho2     ;
  double dT= 0.01 ;

  p     = Mixture_pressure(Temp,density, yi);
  rho2  = MIXTURE_density(Temp+dT,p,yi)     ;

  dhdT  = (MIXTURE_enthalpy(Temp+dT,rho2,P,yi) - MIXTURE_enthalpy(Temp,
           density,P,yi)) /dT ;
    
 
  return dhdT ; /* J/(Kg.K) */
}

double MIXTURE_enthalpy_p(double Temp, double density, double P, double yi[])
{
  double dhdp     ;
  double p        ;
  double rho2     ;
  double dp= 5.0  ;

  p     = Mixture_pressure(Temp,density, yi);
  rho2  = MIXTURE_density(Temp,p+dp,yi)     ;

  dhdp  = (MIXTURE_enthalpy(Temp,rho2,P,yi) - MIXTURE_enthalpy(Temp,density,
           P,yi)) /dp;


  return dhdp ;  /* J/ (Kg.Pascal) */ 
}


/*******************************************************************/
/* Auxiliary Mixture Functions                                      */
/*******************************************************************/


double Mixture_Rgas(double yi[])
{
  double Rgas=0.0 ;
  int i ;

  for (i=0; i<n_specs; ++i)
    Rgas += yi[i]*(UNIVERSAL_GAS_CONSTANT/Mw_i(i)) ;
  
  return Rgas ; 
}

double Mixture_pressure(double Temp, double Rho, double yi[])
{
  double Rgas =  Mixture_Rgas(yi) ;
  
  double P    =  Rho*Rgas*Temp    ;  /* Pressure at Temp & P */

  return P;			/* (Kg/m^3) */
}




/*******************************************************************/
/* Species Property Functions                                      */ 
/*******************************************************************/
double Mw_i(int i)
{
  double mi[20];

  mi[0] = 18.01534 ; /*H2O*/
  mi[1] = 28.01340 ; /*N2 */
  mi[2] = 31.99880 ; /*O2 */
  mi[3] = 44.00995 ; /*CO2*/
  
  return mi[i] ;
}

double Cp_i(double T, double r, int i)
{
  double cpi[20] ;
  
  cpi[0] = 2014.00 ; /*H2O*/
  cpi[1] = 1040.67 ; /*N2 */
  cpi[2] =  919.31 ; /*O2 */
  cpi[3] =  840.37 ; /*CO2*/

  return cpi[i] ;
}

double K_i(double T, double r, int i)
{
  double ki[20] ;
  
  ki[0] = 0.02610 ; /*H2O*/
  ki[1] = 0.02420 ; /*N2 */
  ki[2] = 0.02460 ; /*O2 */
  ki[3] = 0.01450 ; /*CO2*/

  return ki[i] ;
}

double Mu_i(double T, double r, int i)
{
  double mui[20] ;
  
  mui[0] = 1.340E-05 ; /*H2O*/
  mui[1] = 1.663E-05 ; /*N2 */
  mui[2] = 1.919E-05 ; /*O2 */
  mui[3] = 1.370E-05 ; /*CO2*/

  return mui[i] ;
}

double Rgas_i(double T, double r, int i)
{
  double Rgasi ;

  Rgasi =  UNIVERSAL_GAS_CONSTANT/Mw_i(i) ;
  
  return Rgasi ;  
}

double Gm_i(double T, double r, int i)
{
  double gammai ;

  gammai =  Cp_i(T,r,i)/(Cp_i(T,r,i) - Rgas_i(T,r,i));
  
  return gammai ;
  
}



/*******************************************************************/
/* Mixture Functions Structure                                     */ 
/*******************************************************************/
UDF_EXPORT RGAS_Functions RealGasFunctionList = 
{
  MIXTURE_Setup,			/* initialize		  */
  MIXTURE_density,       		/* density	       	  */
  MIXTURE_enthalpy,		        /* enthalpy  		  */
  MIXTURE_entropy,		        /* entropy  		  */
  MIXTURE_specific_heat,		/* specific_heat	  */
  MIXTURE_mw,		        	/* molecular_weight  	  */
  MIXTURE_speed_of_sound,		/* speed_of_sound  	  */
  MIXTURE_viscosity,		        /* viscosity		  */
  MIXTURE_thermal_conductivity,	        /* thermal_conductivity	  */
  MIXTURE_rho_t,			/* drho/dT |const p 	  */
  MIXTURE_rho_p,			/* drho/dp |const T 	  */
  MIXTURE_enthalpy_t,	        	/* dh/dT |const p 	  */
  MIXTURE_enthalpy_p		        /* dh/dp |const T 	  */
};

/*******************************************************************/
/*******************************************************************/



UDRGM Example: Real Gas Model with Volumetric Reactions

This is an example of a UDRGM that has been set up for reacting flow simulations. The example UDF code consists of the following sections:

figure   

In the UDRGM only the mixture species and associated properties are defined. No information about the chemical reactions is required in the UDF. The chemical reaction is set up in a separate step, after the UDF has been compiled and loaded into ANSYS FLUENT. See this section in the separate User's Guide for details. 

/*
 *sccs id: @(#)real_ideal.c	1.10 Copyright 1900/11/09 ANSYS, Inc.
 */

/*
 *      Copyright 1988-1998 ANSYS, Inc.
 *      All Rights Reserved
 *
 *      This is unpublished proprietary source code of ANSYS, Inc.
 *      It is protected by U.S. copyright law as an unpublished work
 *      and is furnished pursuant to a written license agreement.  It
 *      is considered by ANSYS, Inc. to be confidential and may not be
 *      used, copied, or disclosed to others except in accordance with
 *      the terms and conditions of the license agreement.
 */


/*
 *    Warning!!!  Including udf.h is for getting definitions for
 *                  ANSYS FLUENT constructs such as Domain.   You must
 *                  NOT reference any ANSYS FLUENT globals directly from 
 *                  within this module nor link this against any ANSYS 
 *                  FLUENT libs, doing so will cause dependencies on a  
 *                  specific ANSYS FLUENT binary such as fl551.exe and 
 *                  thus won't be version independent.
 */

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

#if RP_DOUBLE
#define SMLL 1.e-20
#else
#define SMLL 1.e-10
#endif

#define NSPECIE_NAME 80 
#define RGASU UNIVERSAL_GAS_CONSTANT /* 8314.34 SI units: J/Kmol/K */
#define PI    3.141592654
/* Here input the number of species in the mixture */
/* THIS IS A USER INPUT */
#define n_specs 5

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

static double ref_p, ref_T;
static char   gas[n_specs][NSPECIE_NAME];

/* static property parameters */
static double cp[5][n_specs]; /* specific heat polynomial coefficients */
static double mw[n_specs];    /* molecular weights */
static double hf[n_specs];    /* formation enthalpy */
static double tcrit[n_specs]; /* critical temperature */
static double pcrit[n_specs]; /* critical pressure */
static double vcrit[n_specs]; /* critical specific volume */
static double nrk[n_specs];   /* exponent n of function a(T) in Redlich-Kwong 
                                 equation of state */
static double omega[n_specs]; /* acentric factor */

/* Static variables associated with Redlich-Kwong Model */
static double rgas[n_specs], a0[n_specs], b0[n_specs], c0[n_specs], 
  bb[n_specs], cp_int_ref[n_specs];

void  Mw();
void  Cp_Parameters(); 
void  Hform();
void  Tcrit();
void  Pcrit();
void  Vcrit();
void  NRK();
void Omega();

double RKEOS_spvol(double temp, double press, int i);
double RKEOS_dvdp(double temp, double density, int i);
double RKEOS_dvdt(double temp, double density, int i);
double RKEOS_H_ideal_gas(double temp, int i);
double RKEOS_specific_heat(double temp, double density, int i);
double RKEOS_enthalpy(double temp, double density, int i);
double RKEOS_entropy(double temp, double density, int i);
double RKEOS_viscosity(double temp, int i);
double RKEOS_thermal_conductivity(double temp, int i);
double RKEOS_vol_specific_heat(double temp, double density, int i);

DEFINE_ON_DEMAND(I_do_nothing)
{
  /*
    This is a dummy function
    must be included to allow for the use of the
    ANSYS FLUENT UDF compilation utility
  */
}

void Mixture_error(int err, char *f, char *msg)
{
  if (err) 
    usersError("Mixture_error (%d) from function: %s\n%s\n",err,f,msg);
}

/*******************************************************************/
/* Mixture Functions                                               */ 
/* These are the only functions called from ANSYS FLUENT Code      */ 
/*******************************************************************/
void MIXTURE_Setup(Domain *domain, cxboolean vapor_phase, char *specielist,
                   int  (*messagefunc)(char *format, ...),
                   void (*errorfunc)(char *format, ...))
{
  /* This function will be called from ANSYS FLUENT after the
     UDF library has been loaded.
     User must enter the number of species in the mixture
     and the name of the individual species.
  */
  
  int i ;
  usersMessage = messagefunc;
  usersError   = errorfunc;
  ref_p        = ABS_P(RP_Get_Real("reference-pressure"),op_pres);
  ref_T        = 298.15 ;

  Message0("\n MIXTURE_Setup: Redlich-Kwong equation of State"
           " with ideal-gas mixing rules \n");
  Message0("\n MIXTURE_Setup: reference-temperature is %f \n", ref_T);

  if (ref_p == 0.0)
    {
      Message0("\n MIXTURE_Setup: reference-pressure was not set by user \n");
      Message0("\n MIXTURE_Setup: setting reference-pressure to 101325 Pa \n");
      ref_p = 101325.0 ;
    }
  /*====================================================*/
  /*=========   User Input Section =====================*/
  /*====================================================*/
  /*
    Define Species name.
    DO NOT use space for naming species
  */
  (void)strcpy(gas[0],"H2O") ; 
  (void)strcpy(gas[1],"CH4") ;
  (void)strcpy(gas[2],"O2")  ; 
  (void)strcpy(gas[3],"CO2") ;
  (void)strcpy(gas[4],"N2")  ; 
  
  /*====================================================*/
  /*=========   End Of User Input Section ==============*/
  /*====================================================*/

  Message0("\n MIXTURE_Setup: RealGas mixture initialization \n");
  Message0("\n MIXTURE_Setup: Number of Species = %d  \n",n_specs);
  for (i=0; i<n_specs; ++i)
    {      
      Message0("\n MIXTURE_Setup: Specie[%d]         = %s  \n",i,gas[i]);
    }

  /*
    concatenate species name into one string
    and send back to fluent
  */
  strcat(specielist,gas[0]);
  for (i=1; i<n_specs; ++i)
    {
      strcat(specielist," ");
      strcat(specielist,gas[i]);
    }

  /* initialize */
  Mw();
  Cp_Parameters(); 
  Hform();
  Tcrit();
  Pcrit();
  Vcrit();
  Omega();
  NRK();

  for (i=0; i<n_specs; ++i)
    {
      rgas[i] = RGASU/mw[i];
      a0[i] = 0.42747*rgas[i]*rgas[i]*tcrit[i]*tcrit[i]/pcrit[i];
      b0[i] = 0.08664*rgas[i]*tcrit[i]/pcrit[i];
      c0[i] = rgas[i]*tcrit[i]/(pcrit[i]+a0[i]/(vcrit[i]*(vcrit[i]+b0[i])))
              +b0[i]-vcrit[i];
      bb[i] = b0[i]-c0[i];
      cp_int_ref[i] = cp[0][i]*log(ref_T)+ref_T*(cp[1][i]+ref_T*(0.5*cp[2][i]
                      +ref_T*(0.333333*cp[3][i]+0.25*cp[4][i]*ref_T))); 
    }
}

double MIXTURE_mw(double yi[])
{
  double MW, sum=0.0 ;
  int i ;

  for (i=0; i<n_specs; ++i)
    sum += yi[i]/mw[i] ;

  MW = 1.0/MAX(sum,SMLL)      ;
  
  return MW;			/* (Kg/Kmol) */
}

double MIXTURE_density(double temp, double P, double yi[])
{
  double den=0.0 ;
  int i ;
    
  for (i=0; i<n_specs; ++i)
    {
      if (yi[i]> SMLL)
   den += yi[i]*RKEOS_spvol(temp, P, i); 
    }

  return 1./den; /* (Kg/m^3) */
}

double MIXTURE_specific_heat(double temp, double density, double P, 
                             double yi[])
{
  double cp=0.0 ;
  int i ;
    
  for (i=0; i<n_specs; ++i)
      if (yi[i]> SMLL)
    cp += yi[i]*RKEOS_specific_heat(temp,mw[i]*density/MIXTURE_mw(yi),i);
  
  return cp;			/* (J/Kg/K) */
}

double MIXTURE_enthalpy(double temp, double density, double P, double yi[])
{
  double h=0.0;
  int i ;

  for (i=0; i<n_specs; ++i)
      if (yi[i]> SMLL)
    h  += yi[i]*RKEOS_enthalpy(temp, mw[i]*density/MIXTURE_mw(yi), i);
  
  return h;			/* (J/Kg) */
}

double MIXTURE_enthalpy_prime(double temp, double density, double P, 
                              double yi[], double hi[])
{
  double h=0.0;
  int i ;

  for (i=0; i<n_specs; ++i)
{
    hi[i] = hf[i]/mw[i] + RKEOS_enthalpy(temp, mw[i]*density/MIXTURE_mw(yi), 
            i);
      if (yi[i]> SMLL)
    h  += yi[i]*(hf[i]/mw[i] + RKEOS_enthalpy(temp, 
          mw[i]*density/MIXTURE_mw(yi), i));
}
  
  return h;			/* (J/Kg) */
}

double MIXTURE_entropy(double temp, double density, double P, double yi[])
{
  double s   = 0.0   ;
  double sum = 0.0 ;
  double xi[n_specs] ;
  int i ;

  for (i=0; i<n_specs; ++i)
    {
      xi[i] = yi[i] / mw[i];
      sum   += xi[i];
    }
  for (i=0; i<n_specs; ++i)
    xi[i] /= sum;

  for (i=0; i<n_specs; ++i)
      if (yi[i]> SMLL)
    s  += yi[i]*RKEOS_entropy(temp,mw[i]*density/MIXTURE_mw(yi), i)-
       UNIVERSAL_GAS_CONSTANT/MIXTURE_mw(yi)* xi[i] * log(xi[i]);

  return s;			/* (J/Kg/K) */
}

double MIXTURE_viscosity(double temp, double density, double P, double yi[])
{
  double mu=0.;
  int i ;

  for (i=0; i<n_specs; ++i)
      if (yi[i]> SMLL)
    mu += yi[i]*RKEOS_viscosity(temp,i);

  return mu;			/* (Kg/m/s) */
}

double MIXTURE_thermal_conductivity(double temp, double density, double P, 
                                    double yi[])
{
  double kt=0.;
  int i ;

  for (i=0; i<n_specs; ++i)
      if (yi[i]> SMLL)
    kt += yi[i]* RKEOS_thermal_conductivity(temp,i);
  
  return kt;			/* W/m/K */
}

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

double MIXTURE_speed_of_sound(double temp, double density, double P, 
                              double yi[])
{
   double dvdp = 0.; 
   double cv = 0.;
   double v = 1./density;   
   int i;
   double cp = MIXTURE_specific_heat(temp, density, P, yi);

   for (i=0; i<n_specs; ++i)
      if (yi[i]> SMLL)
     {
       dvdp += yi[i]*RKEOS_dvdp(temp, mw[i]*density/MIXTURE_mw(yi),i);
       cv += yi[i]*RKEOS_vol_specific_heat(temp, mw[i]*density/MIXTURE_mw(yi), 
             i);
     }

   return sqrt(- cp/cv/dvdp)*v;  
}

/*------------------------------------------------------------*/
/* FUNCTION: MIXTURE_rho_t                                    */
/*------------------------------------------------------------*/

double MIXTURE_rho_t(double temp, double density, double P, double yi[])
{
  double rho_t = 0.;
  int i;

  for (i=0; i<n_specs; ++i)
      if (yi[i]> SMLL)
   rho_t -= yi[i]*density*density*RKEOS_dvdt(temp,
            mw[i]*density/MIXTURE_mw(yi) , i);
  return rho_t;
}

/*------------------------------------------------------------*/
/* FUNCTION: MIXTURE_rho_p                                    */
/*------------------------------------------------------------*/

double MIXTURE_rho_p(double temp, double density, double P, double yi[])
{
  double rho_p = 0.;
  int i;

  for (i=0; i<n_specs; ++i)
      if (yi[i]> SMLL)
   rho_p -= yi[i]*density*density*RKEOS_dvdp(temp, 
            mw[i]*density/MIXTURE_mw(yi), i);
  return rho_p;
}

/*------------------------------------------------------------*/
/* FUNCTION: MIXTURE_enthalpy_t                               */
/*------------------------------------------------------------*/

double MIXTURE_enthalpy_t(double temp, double density, double P, double yi[])
{
   return MIXTURE_specific_heat(temp, density, P, yi);
}

/*------------------------------------------------------------*/
/* FUNCTION: MIXTURE_enthalpy_p                                 */
/*------------------------------------------------------------*/

double MIXTURE_enthalpy_p(double temp, double density, double P, double yi[])
{
   double v = 1./density;
   double dvdt = 0.0;
   int i;

   for (i=0; i<n_specs; ++i)
      if (yi[i]> SMLL)
      dvdt += yi[i]*RKEOS_dvdt(temp,  mw[i]*density/MIXTURE_mw(yi), i);

   return v-temp*dvdt;
}

/*******************************************************************/
/* Species Property Definitions                                    */ 
/*******************************************************************/

void Mw() /* molecular weight */
{  /* Kg/Kgmol */
  mw[0] = 18.01534 ; /*H2O*/
  mw[1] = 16.04303 ; /*CH4*/
  mw[2] = 31.99880 ; /*O2 */
  mw[3] = 44.00995 ; /*CO2*/
  mw[4] = 28.01340 ; /*N2 */
}

void Pcrit() /* critical pressure */
{ /* Pa */
  pcrit[0] = 220.64e5 ; /*H2O*/
  pcrit[1] = 4.48e6 ; /*CH4*/
  pcrit[2] = 5066250.; /*O2 */
  pcrit[3] = 7.3834e6 ; /*CO2*/
  pcrit[4] = 33.98e5 ; /*N2 */
}

void Tcrit() /* critical temperature */
{  /* K */
  tcrit[0] = 647. ; /*H2O*/
  tcrit[1] = 191. ; /*CH4*/
  tcrit[2] = 155.; /*O2 */
  tcrit[3] = 304.; /*CO2*/
  tcrit[4] = 126.2 ; /*N2 */
}

void Vcrit() /* critical specific volume */
{ /* m3/Kg */
  vcrit[0] = 0.003111 ; /*H2O*/
  vcrit[1] = 0.006187 ; /*CH4*/
  vcrit[2] = 0.002294 ; /*O2 */
  vcrit[3] = 0.002136; /*CO2*/
  vcrit[4] = 0.003196; /*N2 */
}

void NRK() /* exponent n of function a(T) in Redlich-Kwong equation of 
              state */
{
  int i;
  for (i=0; i<n_specs; ++i)
    nrk[i]= 0.4986 + 1.1735*omega[i] + 0.475*omega[i]*omega[i];
}
void Omega() /* acentric factor */
{
  omega[0] = 0.348 ; /*H2O*/
  omega[1] = 0.007; /*CH4*/
  omega[2] = 0.021 ; /*O2 */
  omega[3] = 0.225; /*CO2*/
  omega[4] = 0.040; /*N2 */
}

void Hform() /* formation enthalpy */
{
  /*J/Kgmol*/
  hf[0] =  -2.418379e+08; /*H2O*/
  hf[1] =  -74895176. ;   /*CH4*/
  hf[2] = 0. ;            /*O2 */
  hf[3] = -3.9353235e+08 ;/*CO2*/
  hf[4] = 0. ;            /*N2 */
}

void Cp_Parameters( ) /* coefficients of specific heat polynomials */
{  /* J/Kg/K */
  cp[0][0] =  1609.791  ; /*H2O*/
  cp[1][0] = 0.740494; 
  cp[2][0] =-9.129835e-06 ;
  cp[3][0] =-3.813924e-08 ;
  cp[4][0] =4.80227e-12 ;

  cp[0][1] = 872.4671  ; /*CH4*/
  cp[1][1] = 5.305473 ; 
  cp[2][1] = -0.002008295 ;
  cp[3][1] = 3.516646e-07;
  cp[4][1] = -2.33391e-11  ;

  cp[0][2] = 811.1803  ; /*O2 */
  cp[1][2] =0.4108345 ; 
  cp[2][2] =-0.0001750725 ;
  cp[3][2] = 3.757596e-08 ;
  cp[4][2] =-2.973548e-12 ;

  cp[0][3] = 453.577 ; /*CO2*/
  cp[1][3] = 1.65014; 
  cp[2][3] = -1.24814e-3 ;
  cp[3][3] = 3.78201e-7 ;
  cp[4][3] = 0.;

  cp[0][4] = 938.8992 ; /*N2 */
  cp[1][4] = 0.3017911 ; 
  cp[2][4] = -8.109228e-05;
  cp[3][4] =  8.263892e-09  ;
  cp[4][4] = -1.537235e-13 ;
}

/*************************************************************/
/*                                                            */
/* User-Defined Function: Redlich-Kwong Equation of State     */
/*                        for Real Gas Modeling               */
/*                                                            */
/*  Author: Frank Kelecy                                      */
/*    Date: May 2003                                          */
/*Modified: Rana Faltsi                                       */
/*    Date: December 2006                                     */
/*                                                            */
/*                                                            */
/**************************************************************/
/* OPTIONAL REFERENCE (OFFSET) VALUES FOR ENTHALPY AND ENTROPY */

#define H_REF  0.0
#define S_REF  0.0
/*------------------------------------------------------------*/
/* FUNCTION: RKEOS_pressure of species i                      */
/*           Returns pressure given T and density             */
/*------------------------------------------------------------*/

double RKEOS_pressure(double temp, double density, int i)
{
   double v = 1./density;
   double afun = a0[i]*pow(tcrit[i]/temp,nrk[i]);
   return rgas[i]*temp/(v-bb[i])-afun/(v*(v+b0[i]));
}

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

double RKEOS_spvol(double temp, double press, int i)
{
   double a1,a2,a3;
   double vv,vv1,vv2,vv3;
   double qq,qq3,sqq,rr,tt,dd;

   double afun = a0[i]*pow(tcrit[i]/temp,nrk[i]);

   a1 = c0[i]-rgas[i]*temp/press;
   a2 = -(bb[i]*b0[i]+rgas[i]*temp*b0[i]/press-afun/press);
   a3 = -afun*bb[i]/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 = -SIGN(rr)*(pow(sqrt(-dd)+fabs(rr),0.333333));
     vv = (tt+qq/tt)-a1/3.;
   } else {
     if (rr/sqrt(qq3)<-1) {
       tt = PI;
     } else if (rr/sqrt(qq3)>1) {
       tt = 0;
     } 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;
     /*Message0("Three roots %f %f %f \n",vv1, vv2, vv3 );*/
   }

   return vv;
}

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

double RKEOS_dvdp(double temp, double density, int i)
{
   double afun = a0[i]*pow(tcrit[i]/temp,nrk[i]);
   double dterm1,dterm2;
   double v     =  1./ density;

   dterm1   = -rgas[i]*temp*pow((v-b0[i]+c0[i]), -2.0);
   dterm2   = afun*(2.0*v+b0[i])*pow(v*(v+b0[i]),-2.0);

   return 1./( dterm1+dterm2);
}

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

double RKEOS_dvdt(double temp, double density, int i)
{
   double dpdT, dterm1, dterm2;
   double afun = a0[i]*pow(tcrit[i]/temp,nrk[i]);
   double v     = 1./density;

   dterm1   = rgas[i]/(v-b0[i]+c0[i]);
   dterm2   = nrk[i]*afun/((v*(v+b0[i]))*temp);
   dpdT = dterm1+dterm2;

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

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

double RKEOS_Cp_ideal_gas(double temp, int i)
{
  double cpi=(cp[0][i]+temp*(cp[1][i]+temp*(cp[2][i]+temp*(cp[3][i]
             +temp*cp[4][i]))));
  if (cpi<SMLL)
    cpi = 1.0;
  return cpi;
}

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

double RKEOS_H_ideal_gas(double temp, int i)
{
  double h = temp*(cp[0][i]+temp*(0.5*cp[1][i]+temp*(0.333333*cp[2][i]
             +temp*(0.25*cp[3][i]+temp*0.2*cp[4][i]))));
  if (h<SMLL)
    h = 1.0;
  return h;
}

/*-----------------------------------------------------------------*/
/* FUNCTION: RKEOS_vol_specific_heat                               */
/*           Returns constant volume specific heat given T and rho */
/*-----------------------------------------------------------------*/

double RKEOS_vol_specific_heat(double temp, double density, int i)
{
   double afun = a0[i]*pow(tcrit[i]/temp,nrk[i]);
   double v = 1./density;
   double Cv0    = RKEOS_Cp_ideal_gas(temp, i) - rgas[i];
   int np1 = (nrk[i]+1.)/b0[i];
   if (Cv0<SMLL)
     Cv0 = 1.;
   return Cv0 + nrk[i]*np1*afun*log(1.0+b0[i]/v)/temp;

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

double RKEOS_specific_heat(double temp, double density, int i)
{
   double delta_Cp,press,v,dvdt,dadt;
   double afun = a0[i]*pow(tcrit[i]/temp,nrk[i]);

   press = RKEOS_pressure(temp, density, i);
   v = 1./density;
   dvdt = RKEOS_dvdt(temp, density, i);
   dadt = -nrk[i]*afun/temp;
   delta_Cp = press*dvdt-rgas[i]-dadt*(1.+nrk[i])/b0[i]*log((v+b0[i])/v)
            + afun*(1.+nrk[i])*dvdt/(v*(v+b0[i]));

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

}

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

double RKEOS_enthalpy(double temp, double density, int i)
{
   double delta_h ,press, v;

   double afun = a0[i]*pow(tcrit[i]/temp,nrk[i]);

   press = RKEOS_pressure(temp, density, i);
   v = 1./density;
   delta_h = press*v-rgas[i]*temp-afun*(1+nrk[i])/b0[i]*log((v+b0[i])/v);

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

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

double RKEOS_entropy(double temp, double density, int i)
{
   double delta_s,v,v0,dadt,cp_integral;

   double afun = a0[i]*pow(tcrit[i]/temp,nrk[i]);

   cp_integral = cp[0][i]*log(temp)+temp*(cp[1][i]+temp*(0.5*cp[2][i]
                 +temp*(0.333333*cp[3][i]+0.25*cp[4][i]*temp)))
               - cp_int_ref[i];
   if (cp_integral<SMLL)
     cp_integral = 1.0;
   v = 1./density;
   v0 = rgas[i]*temp/ref_p;
   dadt = -nrk[i]*afun/temp;
   delta_s = rgas[i]*log((v-bb[i])/v0)+dadt/b0[i]*log((v+b0[i])/v);

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

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

double RKEOS_viscosity(double temp, int i)
{
   double mu,tr,tc,pcatm;

   tr = temp/tcrit[i];
   tc = tcrit[i];
   pcatm = pcrit[i]/101325.;

   mu = 6.3e-7*sqrt(mw[i])*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,int i)
{
   double cp, mu;

   cp = RKEOS_Cp_ideal_gas(temp, i);
   mu = RKEOS_viscosity(temp, i);

   return (cp+1.25*rgas[i])*mu;
}

/*******************************************************************/
/* Mixture Functions Structure                                     */ 
/*******************************************************************/
UDF_EXPORT RGAS_Functions RealGasFunctionList = 
{
  MIXTURE_Setup,			/* initialize		  */
  MIXTURE_density,       		/* density	       	  */
  MIXTURE_enthalpy,		        /* sensible enthalpy      */
  MIXTURE_entropy,		        /* entropy  		  */
  MIXTURE_specific_heat,		/* specific_heat	  */
  MIXTURE_mw,		        	/* molecular_weight  	  */
  MIXTURE_speed_of_sound,		/* speed_of_sound  	  */
  MIXTURE_viscosity,		        /* viscosity		  */
  MIXTURE_thermal_conductivity,	        /* thermal_conductivity	  */
  MIXTURE_rho_t,			/* drho/dT |const p 	  */
  MIXTURE_rho_p,			/* drho/dp |const T 	  */
  MIXTURE_enthalpy_t,	        	/* dh/dT |const p 	  */
  MIXTURE_enthalpy_p,		        /* dh/dp |const T 	  */
  MIXTURE_enthalpy_prime		/* enthalpy               */
};

next up previous contents index Previous: 8.2.5 User-Defined Scalars
Up: 8.2 Detailed UDF Examples
Next: A. C Programming Basics
Release 12.0 © ANSYS, Inc. 2009-01-14