CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Fluent UDF and Scheme Programming (https://www.cfd-online.com/Forums/fluent-udf/)
-   -   floating point exception (https://www.cfd-online.com/Forums/fluent-udf/202004-floating-point-exception.html)

ctmagic May 17, 2018 22:37

floating point exception
 
Hi, everybody. I tryied to use a non-arrhenius reaction rate in my CFD model.
There are three reactions:
CH4+H2O=3H2+CO
CO+H2O= H2+CO2
CH4+2H2O=4H2+CO2

But I always get the Error reports as below. Could anyone tell me where is the problem?My UDF is also attached.

BTW: How to get the spieces ID? Is the spieces order in material mixture?

Detailed Error Info:
/************************************************** ******/
iter continuity x-velocity y-velocity energy ch4 h2 co co2 time/iter
Stabilizing species-0 to enhance linear solver robustness.

Divergence detected in AMG solver: species-0 Stabilizing species-1 to enhance linear solver robustness.

Divergence detected in AMG solver: species-1 Stabilizing species-2 to enhance linear solver robustness.

Divergence detected in AMG solver: species-2 Stabilizing species-3 to enhance linear solver robustness.

Divergence detected in AMG solver: species-3 Stabilizing temperature to enhance linear solver robustness.

Divergence detected in AMG solver: temperature

WARNING: Invalid cp (0.000000e+00 J/kgK) for mixture-template at temperature -1.#IND00 K


Warning: Mass_Diffusivity: invalid (zero) diffusivity.
Please check material properties.

Error at host: floating point exception

Error at Node 0: floating point exception

Error: floating point exception
Error Object: #f
************************************************** *********/
My UDF
/************************************************/
#include "udf.h"
#include "math.h"
#include "materials.h"
#define CH4 0
#define H2 1
#define CO 2
#define CO2 3
#define H2O 4
#define R_const 8.314
#define cata_density 2355.2
#define cata_porosity 0.51

DEFINE_VR_RATE(reforming_rate,c,t,r,mw,yi,rr,rr_t)
{
double k1, k2, k3, K1, K2, K3,KCH4,KH2,KCO,KH2O;
double Tc = C_T(c,t);
double gas_density = C_R(c,t);
double pCH4,pH2,pCO,pCO2,pH2O;
double DEN;
double r1m,r2m,r3m,r1v,r2v,r3v;
k1 = 9.49e+16*exp(-28879/Tc);
k2 = 4.39e+4*exp(-8074.3/Tc);
k3 = 2.29e+16*exp(-29336.0/Tc);
K1 = 10266.76*exp(-26830 /Tc+30.11);
K2 = exp(4400 /Tc-4.063);
K3 = K1 * K2;
KCH4 = 6.65e-6*exp(4604.28 / Tc);
KH2O = 1.77e+3*exp(-10666.35 / Tc);
KH2 = 6.12e-11*exp(9971.13 / Tc);
KCO = 8.23e-7*exp(8497.71 / Tc);


pCH4=(yi[CH4]*gas_density*R_const*Tc)/(mw[CH4]*1000);
pH2 = (yi[H2]*gas_density*R_const*Tc)/(mw[H2]*1000) ;
pCO = (yi[CO]*gas_density*R_const*Tc)/(mw[CO]*1000) ;
pCO2 = (yi[CO2]*gas_density*R_const*Tc)/(mw[CO2]*1000);
pH2O = (yi[H2O] * gas_density *R_const*Tc) / (mw[H2O]*1000);

DEN = 1 + KCO * pCO + KH2 * pH2 + KCH4 * pCH4 + (KH2O * pH2O / pH2);


r1m = 1000*k1 * (pCH4*pH2O - pow(pH2, 3)*pCO / K1) / (3600*pow(pH2, 2.5)*pow(DEN, 2));
r2m = 1000*k2 * (pCO*pH2O - (pH2*pCO2 / K2)) / (3600*pH2*pow(DEN, 2));
r3m = 1000*k3 * (pCH4*pow(pH2O, 2) - (pow(pH2, 4)*pCO2 / K3)) / (3600*pow(pH2, 3.5)*pow(DEN, 2));
r1v = r1m * cata_density*(1 - cata_porosity);
r2v = r2m * cata_density*(1 - cata_porosity);
r3v = r3m * cata_density*(1 - cata_porosity);

*rr=r1v;

}

obscureed May 21, 2018 08:51

Hi ctmagic,

It is possibly worth paying attention to the error/warning messages about specific heat capacity and diffusivity -- although temperature may have diverged by the time those messages were generated.

I think it is essential to check that your UDF generates correct values. Apart from looking at the results and seeing that they are more or less reasonable, there are no checks apart from the ones you do yourself. So, as a debugging stage, print out every single value in the calculations -- but only for a few cells, or you will be swamped. The cell index, c, is really just an integer, so you can use something like:
Code:

if(c % 100000 == 37) Message("cell=%8d: Tc=%16.8g, DEN=%16.8g\n",c,Tc,DEN);
You might need to vary 100000 according to your model size; 37 is just an arbitrary low integer, and you could vary it to find an interesting or typical cell.

Since you want to print out more or less every variable in the UDF, you could try a "stringification" macro, as in https://gcc.gnu.org/onlinedocs/gcc-4...ification.html. Something like this (but I have not tried it exactly myself):
Code:

#define PRINTVAR(varname) do{ \
  if(c % 100000 == 37) \
    Message("cell=%8d: Tc=%16.8g, " #varname "=%16.8g\n" \
    ,c,Tc,varname); \
  } while(0)
PRINTVAR(k1);
PRINTVAR(k2);
/* etc etc etc */

See also https://gcc.gnu.org/onlinedocs/cpp/S...Semicolon.html

I advise you to cross-check some values in complete detail, by comparison with a spreadsheet or something. You might find, for example, that you have messed up the factor of 1000 when dealing with R/MW[i]. (To be honest, I think you probably have, though I am not sure. Fluent sticks closely to SI units, and the engineering way to accept molecular masses like C=12.0107 etc is to regard them as having units [kg/kmol]. Once you accept [kmol], then UNIVERSAL_GAS_CONSTANT is 8314.34, and it is #define'd as such in materials.h. But it is not safe to twiddle factors of 1000 until the results converge and look OK -- the point is that you need to check the whole thing. It is very easy to type "k1" instead of "K1", for example.) Even after you have checked that you hand over a rate in the correct units -- [kmol/(m3.s)] -- then you might like to check that Fluent is doing its part correctly in a test case.

You might find that the UDF is sometimes called with Tc==0., at least if you run from an initialized start. This is clearly crazy, and Tc is later initialized, but in general it is a good idea to make UDFs robust to crazy inputs -- even just returning all zeros.

By the way, I really like your "#define CH4 0" etc.

I presume you have only one equation, or else you would need to use the clunky tests "if (!strcmp(r->name, "reaction-1"))" etc, as in the help files.

It seems odd that you send back the result in *rr but not *rr_t -- I think this would be correct only for a laminar simulation (which is possible, I suppose).

Good luck!
Ed

ctmagic May 21, 2018 21:20

Hi, obscureed!
Many thanks for your patient reply,it's really helpful! I'll try to debug the UDF according to your suggestion.
Thanks again!

Quote:

Originally Posted by obscureed (Post 693056)
Hi ctmagic,

It is possibly worth paying attention to the error/warning messages about specific heat capacity and diffusivity -- although temperature may have diverged by the time those messages were generated.

I think it is essential to check that your UDF generates correct values. Apart from looking at the results and seeing that they are more or less reasonable, there are no checks apart from the ones you do yourself. So, as a debugging stage, print out every single value in the calculations -- but only for a few cells, or you will be swamped. The cell index, c, is really just an integer, so you can use something like:
Code:

if(c % 100000 == 37) Message("cell=%8d: Tc=%16.8g, DEN=%16.8g\n",c,Tc,DEN);
You might need to vary 100000 according to your model size; 37 is just an arbitrary low integer, and you could vary it to find an interesting or typical cell.

Since you want to print out more or less every variable in the UDF, you could try a "stringification" macro, as in https://gcc.gnu.org/onlinedocs/gcc-4...ification.html. Something like this (but I have not tried it exactly myself):
Code:

#define PRINTVAR(varname) do{ \
  if(c % 100000 == 37) \
    Message("cell=%8d: Tc=%16.8g, " #varname "=%16.8g\n" \
    ,c,Tc,varname); \
  } while(0)
PRINTVAR(k1);
PRINTVAR(k2);
/* etc etc etc */

See also https://gcc.gnu.org/onlinedocs/cpp/S...Semicolon.html

I advise you to cross-check some values in complete detail, by comparison with a spreadsheet or something. You might find, for example, that you have messed up the factor of 1000 when dealing with R/MW[i]. (To be honest, I think you probably have, though I am not sure. Fluent sticks closely to SI units, and the engineering way to accept molecular masses like C=12.0107 etc is to regard them as having units [kg/kmol]. Once you accept [kmol], then UNIVERSAL_GAS_CONSTANT is 8314.34, and it is #define'd as such in materials.h. But it is not safe to twiddle factors of 1000 until the results converge and look OK -- the point is that you need to check the whole thing. It is very easy to type "k1" instead of "K1", for example.) Even after you have checked that you hand over a rate in the correct units -- [kmol/(m3.s)] -- then you might like to check that Fluent is doing its part correctly in a test case.

You might find that the UDF is sometimes called with Tc==0., at least if you run from an initialized start. This is clearly crazy, and Tc is later initialized, but in general it is a good idea to make UDFs robust to crazy inputs -- even just returning all zeros.

By the way, I really like your "#define CH4 0" etc.

I presume you have only one equation, or else you would need to use the clunky tests "if (!strcmp(r->name, "reaction-1"))" etc, as in the help files.

It seems odd that you send back the result in *rr but not *rr_t -- I think this would be correct only for a laminar simulation (which is possible, I suppose).

Good luck!
Ed


azharuddin0613 March 9, 2021 19:52

Hello??

Any update on your issue???

I see similar errors when I compile my udf.

Let me know.

Thanks in advance.

regards,


All times are GMT -4. The time now is 13:16.