[schwablyon]

Hello there,


before I describe my current problem, I'd like to thank you guys for the work that you put into this discussion forum. I've been working on my master thesis for a couple of months now and you've helped me out a couple of times.


But now I've reached the point where I need specific help to a multiphase-UDF reaction rate problem that I couldn't solve so far. I'm trying to model a fluidized bed reactor for methanol synthesis with the use of the Euler-Granular model. Syngas is my primary phase, while the granular phase is designed to be the catalyst material. Currently I'm trying to model the reactions taking place in the reactor. I'm using a kinetic model (by Graaf et. al. from 1988) to determine the reaction rates of the following reactions:


CO2 + 3H2 <-> CH3OH + H2O

CO + 2H2 <-> CH3OH



CO + H2O <-> H2 + CO2



Now, I am trying to implement the DEFINE_VR_RATE that I used for a different reactor. This works so far, but ideally, the reaction rates would only be defined in cells in which catalyst material (solid granular phase) is present. The question is, whether I should use the phase interaction reactions (DEFINE_HET_RXN_RATE) or the volumetric reaction rate (DEFINE_VR_RATE) to solve this problem. I've tried both but I am not getting the results I wished for. (SIGSEV errors due to issues with retrieving the correct phase, instant floating point errors)...


I believe the best way would be to use the volumetric reaction rate and loop over all cells to find out the volume fraction of the granular phase. But doesn't DEFINE_VR_RATE alredy loop over all cells, as it writes the rection rate for all cells? My current code looks like this:

#include"udf.h"



DEFINE_VR_RATE(GraafKinetic, c, t, r, mw, yi, rr, rr_t)
{


#if !PR_HOST



/*initializing variables*/


    int ID = 3;         /* Zone ID for granular phase zone from phases panel */
          Domain *granular_domain = Get_Domain(2);  /* get granular phase domain*/
          
          Thread *thread = Lookup_Thread(granular_domain, ID);  /* get phase thread*/
          real volfr_cat = C_VOF(c, thread);

      if ((volfr_cat > 0.001)){


/*calculating and writing reaction rates*/


}


  #endif  /*PR-HOST*/
}

or

#include"udf.h"



DEFINE_VR_RATE(GraafKinetic, c, t, r, mw, yi, rr, rr_t)
{


#if !PR_HOST



/*initializing variables*/


int ID = 3;         /* Zone ID for granular phase zone from phases panel */
          Domain *mixture_domain = Get_Domain(1);
          Domain *subdomain = DOMAIN_SUB_DOMAIN(mixture_domain,1);
          Thread *thread = Lookup_Thread(subdomain, ID);
        real C_VOF(c,thread) = volfr_cat;

      if (volfr_cat > 0.001)){

/*calculating and writing reaction rates*/


}

  #endif  /*PR-HOST*/
}

I'm getting SIGSEV errors in both cases. Anyone herew ho can help me out?

[souza.emer]

I'm facing exactly the same error with DEFINE_VR_RATE

Have you solved your problem?

Or somebody else could help me?

Here is my UDF:

DEFINE_VR_RATE(homogeneous_reactions, c, t, r, mw, yi, rr, rr_t)
{

double vf_g = C_VOF(c, t);

double rho_g = C_R(c, t);

double mi_g_t = C_MU_T(c, t); /* turbulent viscosity */

double mw_o2 = mw[0];
double mw_co2 = mw[1];
double mw_co = mw[2];
double mw_h2o = mw[3];
double mw_h2 = mw[4];
double mw_volatile_A = mw[5];
double mw_volatile_B = mw[6];
double mw_n2 = mw[7];
double imw_volatile_A = (1. / mw_volatile_A); /* to be used at stoichiometric coefficients calculation */
double imw_volatile_B = (1. / mw_volatile_B);

double yi_o2 = yi[0];
double yi_co2 = yi[1];
double yi_co = yi[2];
double yi_h2o = yi[3];
double yi_h2 = yi[4];
double yi_volatile_A = yi[5];
double yi_volatile_B = yi[6];
double yi_n2 = yi[7];


/* coal A */
double alpha1_A = 2.; /* stoichiometric coefficient chemical equation balance */
double alpha5_A = 0.687; /* stoichiometric coefficient chemical equation balance */

/* coal B */
double alpha1_B = 2.88; /* stoichiometric coefficient chemical equation balance */
double alpha5_B = 1.219; /* stoichiometric coefficient chemical equation balance */

double vrel_d = ((C_U(c, t))*(C_U(c, t)) +
(C_V(c, t))*(C_V(c, t)) +
(C_W(c, t))*(C_W(c, t))); /* Gas velocity magnitude */

double vrel = pow(vrel_d, 0.5); /* REACTS code */

// Turbulent parameters for volatile combustion reaction rates
double kt = pow(0.03 * vrel, 2.);
double et = (pow(kt, 0.5) / 0.01);
double eok = et / kt;

/* Volatile Combustion COAL A */
double volrr_term1_A = 23.6*pow((mi_g_t * eok / rho_g), 0.25);
double volrr_termvol_A = C_DPMS_CONCENTRATION(c, t) * yi_volatile_A;
double volrr1_termo2_A = (vf_g * yi_o2 * mw_volatile_A) / (alpha1_A*mw_o2);
double volrr2_termco2_A = (vf_g * yi_co2 * mw_volatile_A) / (alpha5_A*mw_co2);

/* Volatile Combustion COAL B */
double volrr_term1_B = 23.6*pow((mi_g_t * eok / rho_g), 0.25);
double volrr_termvol_B = C_DPMS_CONCENTRATION(c, t) * yi_volatile_B;
double volrr1_termo2_B = (vf_g * yi_o2 * mw_volatile_B) / (alpha1_A*mw_o2);
double volrr2_termco2_B = (vf_g * yi_co2 * mw_volatile_B) / (alpha5_A*mw_co2);

if (!strcmp(r->name, "reaction-5"))
{
/* Reaction 5 - Volatiles + alpha1 O2 = alpha2 CO2 + alpha3 H2O + alpha4 N2 */

double volrr1_A = volrr_term1_A * eok * MIN(volrr_termvol_A, volrr1_termo2_A);

*rr = volrr1_A;
/* *rr_t = *rr; */

}


else if (!strcmp(r->name, "reaction-6"))
{
/* Reaction 4 - Volatiles + alpha5 CO2 = alpha6 CO + alpha7 H2 + alpha8 N2 */

/* Nogami alphas - alpha5 = 0.497, alpha6 = 1.016, alpha7 = 0.425, alpha8 = 0.005 */

double volrr2_A = volrr_term1_A * eok * MIN(volrr_termvol_A, volrr2_termco2_A);

*rr = volrr2_A;
/* *rr_t = *rr; */

}

else if (!strcmp(r->name, "reaction-7"))
{
/* Reaction 5 - Volatiles + alpha1 O2 = alpha2 CO2 + alpha3 H2O + alpha4 N2 */

double volrr1_B = volrr_term1_B * eok * MIN(volrr_termvol_B, volrr1_termo2_B);

*rr = volrr1_B;
/* *rr_t = *rr; */

}


else if (!strcmp(r->name, "reaction-8"))
{
/* Reaction 4 - Volatiles + alpha5 CO2 = alpha6 CO + alpha7 H2 + alpha8 N2 */

/* Nogami alphas - alpha5 = 0.497, alpha6 = 1.016, alpha7 = 0.425, alpha8 = 0.005 */

double volrr2_B = volrr_term1_B * eok * MIN(volrr_termvol_B, volrr2_termco2_B);

*rr = volrr2_B;
/* *rr_t = *rr; */

}

}