DMFC UDF problems
 

I'm modeling DMFC and the reaction source terms is modeled by UDFs. When I run the simulation, it prompts "Error: Divergence detected in AMG solver: x-momentum", I reduce the relaxation factor, it has the same problem, but if I remove the udf, it can work. So I think maybe some problems arise in the UDFs. I use Message function to output the the result of calculation (italic part of the code), some results are "-1.#IND0000,0.99800210,-1.#IND0000,0.99589283", I don't know what's happened. Can anyone help me? Some codes are attached:

#include "udf.h"

#define R UNIVERSAL_GAS_CONSTANT /*gas constant 8314.3 [J/kmol.K] this is known within Fluent*/
#define F 9.65e+7 /* Faraday constant C/kmol */
#define E_0 0.7 /*cell potential [V]*/
#define E_cell 0.6 /*output potential [V]*/

#define mw_ch 32.0 /*molecular weight of methanol [kg/kmol] */
#define mw_co2 44.0 /*molecular weight of carbon dioxide [kg/kmol] */
#define mw_h2o 18.0 /*molecular weight of oxygen [kg/kmol]*/
#define mw_o2 32.0 /*molecular weight of water [kg/kmol] */


#define t_cl 0.02e-3 /* Cathalyst layer thickness [m] */
#define t_mem 0.183e-3 /*membrane thickness [m]*/

#define nd_w 2.5 /*electro-osmotic drag coefficient for water*/
/******The followings are coefficients(need to be adjusted)******/
#define C1 1.723e-7 /*R/(alpha_a*F)*/
#define C2 0.011 /*C_ch_ref/i_a_ref*/
#define C3 8.616e-8 /*R/(alpha_c*F)*/
#define C4 0.011 /*C_O2_ref/i_c_ref*/
#define A_CV 35.65 /*anode catalyst layer surface to volume ratio [1/m]*/
/******The followings are coefficients(need to be adjusted)******/

/*Anode local current density solved by Secant method iteratively */
real Local_current_density_anode(cell_t c, Thread *t)
{
real x0,x1,x2,err;
x0=100.0; /*initial value X0*/
x1=200.0; /*initial value X1*/
err = fabs(x1-x0); /*error*/
while(err>0.000001) /*0.000001 is the convergency condition*/
{
x2= x1-F_anode_local_current_density(c,t,x1)*(x1-x0)/(F_anode_local_current_density(c,t,x1)-F_anode_local_current_density(c,t,x0));//Secand method
x0=x1;
x1=x2;
err=fabs(x1-x0);
if(err<=0.000001)
{
break;
}
}
return x2;
}

/*function of anode local current density*/
real F_anode_local_current_density(cell_t c, Thread *t, real x)
{
real annodeOverpotential, anodeReferenceCurrentDensity,cathodeOverpotential, cathodeReferenceCurrentDensity,anodeLocalCurrentDe nsity;

anodeReferenceCurrentDensity = 1.0*exp((35570.0/R)*(1.0/353.0-1.0/C_T(c,t))); /*anode reference current density*/
annodeOverpotential = C1*C_T(c,t)*log(C2*x/(anodeReferenceCurrentDensity*C_T(c,t)*C_R(c,t)/mw_ch)); /*anode overpotential*/

cathodeReferenceCurrentDensity = 1.0*exp((73200.0/R)*(1.0/353.0-1.0/C_T(c,t))); /*cathode reference current density*/
cathodeOverpotential = C3*C_T(c,t)*log(C4*((1.0+6.0*nd_w*(C_YI(c,t,0)*(C_ YI(c,t,0)*mw_ch+C_YI(c,t,1)*mw_co2+C_YI(c,t,2)*mw_ h2o)/mw_ch))*x+6.0*F*4.012e-

13*exp(0.024312*C_T(c,t))*pow(0.2,1.5)*C_R(c,t)*C_ YI(c,t,0)/(t_mem*mw_ch))/(cathodeReferenceCurrentDensity*0.21*1.29/mw_o2));/*cathode overpotential; 1.29 is air density [kg/m^3]*/

anodeLocalCurrentDensity = 7.3*exp(1268.0*(1.0/298.0-1.0/C_T(c,t)))*(E_0-E_cell-annodeOverpotential-cathodeOverpotential)/t_mem-x;
Message("%5.12f,%5.8f,%5.12f,%5.8f\n",annodeOverpo tential,anodeReferenceCurrentDensity,cathodeOverpo tential, cathodeReferenceCurrentDensity);
return anodeLocalCurrentDensity;
}

/*methanol source term1, caused by the reaction*/
DEFINE_SOURCE(meoh_source1,c,t,dS,eqn)
{
real source11;
source11 = 0.0+(-mw_ch*Local_current_density_anode(c,t)*A_CV/(6.0*F));
dS[eqn] = 0.0;
return source11;
}

Can anyone help me to solve these problems? Please!