[longbk]

Dear all,

I have been conducting the project related to modeling the hydrogen storage based on metal hydride. In my case, hydrogen is injected into the reactor, then hydrogen will react with metal (porous material) to form metal hydride. The amount of absorbed hydrogen is represented by "sink term" which is covered in my UDF.
Here is my UDF code:

#include"udf.h"

/* Input Boundary Conditions */
#define T_i 293.0 /* Inlet temp in K */
#define T_f 293.0 /* Coolant fluid temp in Kelvin [K] */
#define P_i 500000.0 /* Inlet pressure in pascals */
#define x_i 0.0 /* Initial value of H/M for reaction to START */
#define x_f 1.0 /* Final value of H/M for reaction to STOP @ saturation */

/* Properties of Metal Alloy [LaNi5] */
#define rho_m 8394.0 /* Density in kg/m^3 */
#define Cp_m 419.0 /* Sp heat in J/kg-K */
#define K_m 1.6 /* Thermal conductivity in W/m-K */
#define por 0.5 /* Porosity */
#define NA 6.0 /* Number of atoms in Mischmetal alloy [MmNi4.6Al0.4] */
#define M_m 434.0 /* molecular weight of mischmetal alloy [MmNi4.6Al0.4] in g/mol; [LaNi5 = 434] */
#define rho_ss 8394.0 /* Density of MH @ saturation kg/m^3 */
#define rho_s 8280.0 /* Density of MH kg/m^3 */
#define E_a 21179.6 /* Activation energy in J/mol H2 */
#define per 1.0e-8 /* Permeability */
#define DELTA_S 107.2 /* Entropy of Formation J/mol H2-K */
#define DELTA_H 30780.0 /* Enthalpy of Formation J/mol H2 */
#define wt 0.013 /* Hydrogen weight Storage Capacity 1.3% of Mischmetal */
#define A_vantHoff 17.738
#define B_vantHoff 3704.6

/* Properties of Hydrogen */
#define K_g 0.127 /* Thermal Conductivity in W/m-K */
#define Cp_g 14283 /* Sp heat in J/kg-K */
#define rho_g 0.0838 /* Density in kg/m^3 */
#define M_g 2.016 /* Molecular weight in g/mol */

/* CONSTANTS */
#define R_u 8.314 /* Universal gas Constant J/mol-K */
#define k_a 59.187 /* Reaction Constant for absorption in (sec)^-1 */

real x_t_temporary;
real rate_temporary;
real heat_rate_temporary;
//////////////////////////////////////////////////////////////////////////////////////////////////
DEFINE_INIT(my_init_fuc,d)
{
cell_t c;
Thread *t;
thread_loop_c(t,d)
{
begin_c_loop_all(c,t)
{
C_UDMI(c,t,0) = x_i;
C_UDMI(c,t,1) = 0.0;
C_UDMI(c,t,2) = 0.0;
}
end_c_loop_all(c,t)
}
}

DEFINE_SOURCE(mass_source_edited1,c,t,dS,eqn)
{
real tp;
real rate;
real P_eq;
real cond;
int curr_N_timestep;
real curr_timestep_size;
real curr_rho;
real x_t;
real dxdt;
real P_g = C_P(c,t);
real vol = C_VOLUME(c,t);
tp = C_T(c,t);
P_eq = pow(10,5)*exp((DELTA_S/R_u)-(DELTA_H/(R_u*tp)));
cond = P_g/P_eq;
curr_N_timestep = N_TIME; //Number of timestep
curr_timestep_size = CURRENT_TIMESTEP; //current timestep size

if(cond > 1.0)
{
x_t = 1.0 - exp((-1)*(rho_ss-rho_s)*(1-por)*k_a*exp(-E_a/(R_u*tp))*((P_g/P_eq)-1)*CURRENT_TIME);
dxdt = (x_t - C_UDMI(c,t,0))/curr_timestep_size;
if(CURRENT_TIME > 0)
{
rate = dxdt*(rho_ss-rho_s)*(1-por);
}
dS[eqn] = 0.0;
}
else
{
rate = 0.0;
dS[eqn] = 0.0;
}

x_t_temporary = x_t;

rate_temporary = rate;

return rate;
}

DEFINE_SOURCE(heat_source_edited1,c,t,dS,eqn)
{
real heat_rate;
real tp;
real rate;
real P_eq;
real cond;
int curr_N_timestep;
real curr_timestep_size;
real x_t;
real dxdt;
real P_g = C_P(c,t);
real vol = C_VOLUME(c,t);
tp = C_T(c,t);

P_eq = pow(10,5)*exp((DELTA_S/R_u)-(DELTA_H/(R_u*tp)));
cond = P_g/P_eq;
curr_N_timestep = N_TIME;
curr_timestep_size = CURRENT_TIMESTEP;

if(cond > 1.0)
{
x_t = 1.0 - exp((-1)*(rho_ss-rho_s)*(1-por)*k_a*exp(-E_a/(R_u*tp))*((P_g/P_eq)-1)*CURRENT_TIME);
dxdt = (x_t - C_UDMI(c,t,0))/curr_timestep_size;
if(CURRENT_TIME > 0)
{
rate = dxdt*(rho_ss-rho_s)*(1-por);
heat_rate = (1000.0*DELTA_H/M_g)*(rate);
}

dS[eqn] = 0.0;

}
else
{
heat_rate = 0.0;
dS[eqn] = 0.0;
}

heat_rate_temporary = heat_rate;
return heat_rate;
}

DEFINE_EXECUTE_AT_END(ThisRunsAtEndOfTimestep)
{

Domain *d;
Thread *t;
cell_t c;
d = Get_Domain(1);

thread_loop_c(t,d)
{
if(FLUID_THREAD_P(t))
{
begin_c_loop(c,t)
{
C_UDMI(c,t,0) = x_t_temporary;
C_UDMI(c,t,1) = rate_temporary;
C_UDMI(c,t,2) = heat_rate_temporary;
}
end_c_loop(c,t)

}
}
}

I tried to run this UDF many times but it seems not working... Anyone could help me?

Regards,