Learning Forum

[lai840713]

#include "udf.h"

#include "math.h"

#include "unsteady.h"

/* Material Properties */

 

#define R_p_CO2 3.75E-4 /*m*/

#define porosity 0.46

#define roh_ads_CO2 1138 /*kg/m^3*/

#define w_m_CO2 0.4423 /*kg/kg*/

#define K0_CO2 6.91E-5 /*mbar^-1*/

#define del_h_ads_CO2 -15000 /*J/kg*/

#define R_u 8.314 /*J/mol/K*/

#define MW_ads 0.01 /*kg/mol*/

#define n_CO2 2

#define D_s0_CO2 2.54E-4 /*m^2/s*/

#define E_a_CO2 42000.0 /*J/mol*/

 

DEFINE_UDS_UNSTEADY(my_uds_unsteady,c,t,i,apu,su)

{

real physical_dt, vol, rho, phi_old;

physical_dt = RP_Get_Real("physical-time-step");

vol = C_VOLUME(c,t);

rho = 1.0; /*C_R_M1(c,t);*/

*apu = -rho*vol / physical_dt; /*implicit part*/

phi_old = C_STORAGE_R(c,t,SV_UDSI_M1(i));

*su = rho*vol*phi_old/physical_dt; /*explicit part*/

}

 

DEFINE_SOURCE(adsorption_CO2,c,t,dS,eqn)

{

real p_op,P_abs,const1,source,w_eq;

real k1,k2,k_m;

p_op = RP_Get_Real("operating-pressure");

P_abs = (C_P(c,t) + p_op) * 1.0E-2; /*mbar*/

const1 = K0_CO2 * exp(del_h_ads_CO2 * MW_ads / (R_u * C_T(c,t))) * P_abs;

w_eq = w_m_CO2 * const1 / pow(1.0 + pow(const1,n_CO2), 1.0/n_CO2);

k1 = 15.0 * D_s0_CO2 / pow(R_p_CO2,2.0); /*1/s*/

k2 = E_a_CO2 / R_u; /*K*/

k_m = k1 * exp(-k2 / C_T(c,t)); /*1/s*/

source = (1.0 - porosity) * k_m * (w_eq - C_UDSI(c,t,0));

dS[eqn] = (-k_m) * (1.0 - porosity);

return source;

}

 

 

DEFINE_SOURCE(continuity_source_CO2,c,t,dS,eqn)

{

real p_op,w_eq,P_abs,const1,source;

real k1,k2,k_m;

p_op = RP_Get_Real("operating-pressure");

P_abs = (C_P(c,t) + p_op) * 1.0E-2;/*mbar*/

const1 = K0_CO2 * exp (del_h_ads_CO2 * MW_ads / (R_u * C_T(c,t))) * P_abs;

w_eq = w_m_CO2 * const1 / pow(1.0 + pow(const1,n_CO2), 1.0/n_CO2);

k1 = 15.0 * D_s0_CO2 / pow(R_p_CO2,2.0); /*1/s*/

k2 = E_a_CO2 / R_u; /*K*/

k_m = k1 * exp(-k2 / C_T(c,t)); /*1/s*/

source = (1.0 - porosity)* roh_ads_CO2 * k_m * (w_eq - C_UDSI(c,t,0));

dS[eqn] = -(1.0 - porosity)* roh_ads_CO2 * k_m;

return source;

}

 

DEFINE_SOURCE(energy_source_CO2,c,t,dS,eqn)

{

real p_op,w_eq,P_abs,const1,source;

real k1,k2,k_m;

p_op = RP_Get_Real("operating-pressure");

P_abs = (C_P(c,t) + p_op) * 1.0E-2;/*mbar*/

const1 = K0_CO2 * exp (del_h_ads_CO2 * MW_ads / (R_u * C_T(c,t))) * P_abs;

w_eq = w_m_CO2 * const1 / pow(1.0 + pow((const1),n_CO2), 1.0/n_CO2);

k1 = 15.0 * D_s0_CO2 / pow(R_p_CO2,2.0); /*1/s*/

k2 = E_a_CO2 / R_u; /*K*/

k_m = k1 * exp(-k2 / C_T(c,t)); /*1/s*/

source = -25 *(1.0 - porosity)* roh_ads_CO2 * del_h_ads_CO2 * k_m * (w_eq - C_UDSI(c,t,0));

dS[eqn] = 25 *(1.0 - porosity)* roh_ads_CO2 * del_h_ads_CO2 * k_m;

return source;

}

 

DEFINE_SOURCE(energy_sourcesecond_CO2,c,t,dS,eqn)

{

real p_op,P_abs,const1,source,w_eq;

real k1,k2,k_m;

p_op = RP_Get_Real("operating-pressure");

P_abs = (C_P(c,t) + p_op) * 1.0E-2; /*mbar*/

const1 = K0_CO2 * exp(del_h_ads_CO2 * MW_ads / (R_u * C_T(c,t))) * P_abs;

w_eq = w_m_CO2 * const1 / pow(1.0 + pow(const1,n_CO2), 1.0/n_CO2);

k1 = 15.0 * D_s0_CO2 / pow(R_p_CO2,2.0); /*1/s*/

k2 = E_a_CO2 / R_u; /*K*/

k_m = k1 * exp(-k2 / C_T(c,t)); /*1/s*/

source = 100000 * (w_eq - C_UDSI(c,t,0)) * (w_eq - C_UDSI(c,t,0));

dS[eqn] = -200000 * (w_eq - C_UDSI(c,t,0));

return source;

}

[tsiriaks]

I would recommend to post this question under Physics Simulation  -> Fluid Dynamics  section.

You might get a proper attention there.

Thanks,

Win

[Rob]

I'd also suggest posting what it does: staff are not going to debug the code, but may be able to give pointers based on what's (not) happening.