[David Artemyev]

Hello, everyone!

I am currently working on the CFD modelling of the process steam methane reforming. I wrote a user-defined function which contains the kinetics of SMR according to the following article: https://www.sci-hub.ru/10.1016/j.cej...9d62o743096007. The UDF file contains the following code:

#include "udf.h"
DEFINE_VR_RATE(msr, c, t, r, mw, yi, rr, rr_t)
{
double y_co2 = yi[1];
double y_n2 = yi[5];
double y_co = yi[2];
double y_h2 = yi[4];
double y_ch4 = yi[0];
double y_h2o = yi[3];
double Nsum, R, TP, temp, k1, k2, k3, Keq1, Keq2, Keq3, K_h2, K_co, K_h2o, K_ch4, DEN, r1, r2, r3, dcb;
/*calculate species i in the unit of kgmol i/kg mix*/
y_co2 *= 1 / mw[1];
y_n2 *= 1 / mw[5];
y_co *= 1 / mw[2];
y_h2 *= 1 / mw[4];
y_ch4 *= 1 / mw[0];
y_h2o *= 1 / mw[3];
/*total mole number per kg mix */
Nsum = y_ch4 + y_h2o + y_h2 + y_co + y_co2 + y_n2;
/*calculate mole fraction of species i in the unit of kgmol
i/kgmol mix*/
y_ch4 *= 1 / Nsum;
y_h2o *= 1 / Nsum;
y_h2 *= 1 / Nsum;
y_co *= 1 / Nsum;
y_co2 *= 1 / Nsum;
y_n2 *= 1 / Nsum;
/*gas constant, J/molK */
R = 8.314;
/*equilibrium constant calculated by the empirical formula*/
temp = C_T(c, t);
Keq1 = 10266.76*exp(-26.83 / temp + 30.11)*1e-4; /*bar^2*/
Keq2 = exp(4400 / temp - 4.063);
Keq3 = Keq1 * Keq2; /*bar^2*/
/*reaction rate constant*/
k1 = 4.225e+15*exp(-240100.0 / (R*temp)); /*kmol*bar^-0.5*kg^-1*h^-1*/
k2 = 1.955e+6*exp(-67130.0 / (R*temp));
k3 = 1.02e+15*exp(-243900.0 / (R*temp));
/*adsorption coefficients*/
K_h2 = 6.12e-9*exp(82900.0 / (R*temp));
K_co = 8.23e-5*exp(70650.0 / (R*temp));
K_h2o = 1.77e+5*exp(-88680.0 / (R*temp));
K_ch4 = 6.65e-4*exp(38280.0 / (R*temp));
/*total pressure in the cell near wall, bar*/
TP = C_P(c, t) / (1.0e+5);
/*constant DEN factor for rate equations*/
DEN = 1.0 + TP * (K_ch4*y_ch4 + K_co * y_co + K_h2 * y_h2) + K_h2o * y_h2o / y_h2;
/*unit conversion to kmol/(m3.s) with density catalyst bulk*/
dcb = 0.305556;
/*reaction rate equations*/
r1 = dcb * ((k1 / (pow(TP*y_h2, 2.5)))*((pow(TP, 2.0)*y_ch4*y_h2o) - ((pow(TP, 4.0)*pow(y_h2, 3.0)*y_co) / Keq1))) / pow(DEN, 2.0);
r2 = dcb * ((k2 / (TP*y_h2))*((pow(TP, 2.0)*y_co*y_h2o) - ((pow(TP, 2.0))*y_h2*y_co2) / Keq2)) / pow(DEN, 2.0);
r3 = dcb * ((k3 / (pow(TP*y_h2, 3.5)))*((pow(TP, 3.0)*y_ch4*pow(y_h2o, 2.0)) - ((pow(TP, 5.0)*pow(y_h2, 4.0)*y_co2) / Keq3))) / pow(DEN, 2.0);

if (!strcmp(r->name, "r1"))
{
*rr = r1;
}
else if (!strcmp(r->name, "r2"))
{
*rr = r2;
}
else if (!strcmp(r->name, "r3"))
{
*rr = r3;
}

}

However, after hybrid initialization I get the floating point exception error. Can I count on you for your help in solving the problem?

[DreyCal]

Hello.

I'm unsure if your problem has been resolved since it has been a long while, but I would like to reply regardless so it could possibly help others in the same predicament.

There are many possible causes of the floating point exception error, and diagnosing its cause is not always straightforward. This error can be caused but not limited to by any of the following:

· Problematic mesh
· Solution setup
· The UDF code

For starters, I highly suggest you use the “Check Case” option, which can be found by going to the “Solution” tab, then find it under the “Run Calculation” pane (for ANSYS 2023). The “Check Case” option is very helpful since it check for problematic inputs acroos your whole solution from the top down from the mesh to your solution setup. In my case, I also encountered the floating point exception error, and it was able to diagnose that I was missing an axis of symmetry in the x-axis. The issue was I set-up my 2d simulation for methanol steam reforming for symmetry with respect to the y-axis, and I further read that axisymmetric cases in ANSYS fluent only recognize the axis of symmetry along the x-axis. My solution was simply to rotate the whole domain by 90 degrees so my axis of symmetry was now oriented at the x-axis instead of the y-axis.

If the problem is being caused by the problematic mesh, then you should check mesh quality metrics, such as element quality, orthogonality, skewness, and aspect ratio. There may have been some error when your tried to mesh your domain, and some elements may have been defeatured as a result which could cause the floating point exception error. For a meaning for mesh quality metrics, you should check this link out:

https://featips.com/2022/11/21/ansys...ics-explained/

Then just enable adaptive sizing for your mesh, then mess with the different options bit by bit until you get good values for mesh quality. But seriously, the easiest way to judge if your mesh is good is as long as it looks like a square or an equilateral/isosceles-right triangle for 2d mesh and a cube or regular polyhedron for 3d mesh. However, improving mesh quality more commonly resolves issues of solution convergence than floating point exception error.

If the problem is being caused by your solution setup, well then your best friends are the warning messages and error messages in the console and the “Check Case” option that I mentioned. Check the messages from the console and modify any related parts of your solution, and check any recommendations given by the “Check Case” option. If the “Check Case” option gives something like “General: check the mesh quality” despite having good quality mesh and you’re still encountering the problem, then you’ll just have to brace yourself and do the ole’ trial and error and scanning your whole setup for the problem. Check your species transport setup, cell zone conditions, boundary conditions, and other pertinent setups. Sadly, there are no shortcuts.

If the problem is being caused by your UDF code, then you can diagnose by assigning constant values to any variables in your UDF like replacing the temperature function “C_T(c,t)” with a constant value of 1,000. The floating point exception error usually comes up when divergence is detected in your soluton, and this means something and somewhere in your solution is approaching infinity. You could check each of your variables by replacing them with a constant. If the floating point exception error ceases after this substitution, then you’ll have to revise the equation where the problematic variables were modified with your assigned constants. You could also just do the ole check your whole equation system and identify the problem yourself. In my case, I noticed that I forgot to sum-up my mole fractions, which was the cause of the floating point exception error.

Additionally, the issue in your UDF could be more about the strcmp command that compares the strings in your UDF to the strings in ANSYS. What strcmp in the UDF does it compare the strings of the reaction names in the UDF and in the ANSYS setup, and this means that the reaction names in your UDF and in the ANSYS setup need to have the EXACT SAME NAME. In the case for this post, the strings of the reaction names are “r1”, “r2”, and “r3”, while the default reaction names in the inside the species transport setup of ANSYS are “reaction-1”, “reaction-2”, and “reaction-3”. If these names in the UDF and the ANSYS setup do not match, then ANSYS will return a false value and load the default equations and inputs in your species transport setup related to your reaction rate and completely ignore the reaction rate equations in your UDF, causing the floating point exception error. For simplicity, just name your reaction names in the UDF as “reaction-#” and ensure that the total number of reactions are correct in the ANSYS setup (alongside the correct stoichiometric coefficients).

I hope this helps others who encountered the same issues with steam reforming.