[Mohammadreza01]

Hello everyone,
I am trying to simulate a CO2 separation process from flue gas by a membrane. But I faced some problems in writing the UDF. I have spent two weeks on this problem, but couldn't solve it. I would greatly appreciate any help.
My geometry consists of two rectangular channels, a shell and a tube, and the problem is axisymmetric (please see the attached figure). The flue gas (CO2, N2, and O2) enters through the velocity inlet and exits from the pressure outlet on the shell side. A sweep gas consisting of CO2 (with traces of O2) enters through the velocity inlet and exits from the pressure outlet on the tube side. I want to calculate the CO2 permeation flux based on the temperature and partial pressures of CO2 and O2 in the shell and tube sides. I found a relation in the literature for doing this, and I wrote the following UDF. Also, I specified the velocity inlet boundary condition for the bottom boundary of the shell side and the top boundary of the tube side, where CO2 permeation occurs. The problem with my UDF is that the compiler does not go into the second begin_f_loop. I can't figure out why this is happening. Does anybody have any idea?



#include "udf.h"


// GEOMETRY CONSTANTS
#define L 0.015 // length of the domain [m]
#define N 170 // number of mesh cells along the domain

// BOUNDARY ID
#define IDmss 13 // boundary ID for membrane_shell_side
#define IDmts 12 // boundary ID for membrane_tube_side

// SPECIES ID
#define idCO2 0
#define idO2 1
#define idN2 2

// MODEL CONSTANTS
#define Po 101325 // operating pressure [Pa]
#define e_t_MC 0.145
#define e_t_SO 0.48
#define Zv 2 // oxygen vacancy charge
#define Zc -2 // carbonate ion charge
#define Zh 1 // electron hole charge
#define F 9.649E4 // Faraday's constant [C/mol]
#define R 8.3145 // ideal gas constant [J/mol-K]
#define Mc 0.02004 // molecular weight of the carbonate mixture [kg/mol]
#define Av 3.795E8 // [K/Ohm-m]
#define Ev 102415 // [J/mol]

// OTHER CONSTANTS
#define PI 3.14159265358979
#define Mco2 0.04400995 // molecular weight of CO2 [kg/mol]
#define Mo2 0.0319988 // molecular weight of O2 [kg/mol]
#define Mn2 0.0280134 // molecular weight of N2 [kg/mol]

// VARIABLES
static real Jco2[N]; // CO2 flux through membrane
static real Jo2[N]; // O2 flux through membrane
static real XmCO2; // mass fraction of CO2 in permeate
static real XmO2; // mass fraction of O2 in permeate
static real Rt; // tube radius [m]
static real Rm; // membrane radius [m]
static real hs; // boundary cell height at the shell side
static real ht; // boundary cell height at the tube side
static real Ts; // temperature at shell side
static real PsCO2; // partial pressure of CO2 at shell side
static real PsO2; // partial pressure of O2 at shell side
static real Tt; // temperature at tube side
static real PtCO2; // partial pressure of CO2 at tube side
static real PtO2; // partial pressure of O2 at tube side


DEFINE_ADJUST(flux_calculator, domain)
{
Thread *ts,*tt,*t0,*t1;
face_t f,ff;
cell_t c,c0,c1;
domain = Get_Domain(1);
ts = Lookup_Thread(domain,IDmss);
tt = Lookup_Thread(domain,IDmts);

real x0[ND_ND],xc0[ND_ND],x1[ND_ND],xc1[ND_ND];
real XsCO2,XsO2,XsN2,YsCO2,YsO2,YsN2;
real XtCO2,XtO2,XtN2,YtCO2,YtO2,YtN2;
real Cc; // carbonate ion concentration [mol/m^3]
real rho_c; // density of the carbonate mixture [kg/m^3]
real Dc; // carbonate ion diffusivity [m^2/s]
real Sv; // oxygen vacancy conductivity [1/Ohm-m]
real Sh; // electron hole conductivity [1/Ohm-m]
real H; // membrane thickness
real term1,term2;


// GETTING DATA FROM SHELL SIDE
FILE *fp;
fp = fopen ("Output01.txt","a");

int ii = 0;
begin_f_loop(f,ts)
{
c0 = F_C0(f,ts);
t0 = THREAD_T0(ts);
C_CENTROID(xc0,c0,t0);
F_CENTROID(x0,f,ts);

// mass fractions at the shell side
XsCO2 = C_YI(c0,t0,idCO2);
XsO2 = C_YI(c0,t0,idO2);
XsN2 = C_YI(c0,t0,idN2);

// calculating mole fractions at the shell side
YsCO2 = XsCO2/Mco2/(XsCO2/Mco2 + XsO2/Mo2 + XsN2/Mn2);
YsO2 = XsO2/Mo2/(XsCO2/Mco2 + XsO2/Mo2 + XsN2/Mn2);
YsN2 = XsN2/Mn2/(XsCO2/Mco2 + XsO2/Mo2 + XsN2/Mn2);

Ts = F_T(f,ts);
PsCO2 = (F_P(f,ts) + Po)*YsCO2;
PsO2 = (F_P(f,ts) + Po)*YsO2;

Rm = x0[1];
hs = (xc0[1]-x0[1])*2;

// GETTING DATA FROM THE TUBE SIDE
begin_f_loop(ff,tt)
{
F_CENTROID(x1,ff,tt);
if (fabs(x0[0]-x1[0])/x0[0] < 1E-3)
{
c1 = F_C0(ff,tt);
t1 = THREAD_T0(tt);
C_CENTROID(xc1,c1,t1);

// mass fractions at the tube side
XtCO2 = C_YI(c1,t1,idCO2);
XtO2 = C_YI(c1,t1,idO2);
XtN2 = C_YI(c1,t1,idN2);

// calculating mole fractions at the tube side
YtCO2 = XtCO2/Mco2/(XtCO2/Mco2 + XtO2/Mo2 + XtN2/Mn2);
YtO2 = XtO2/Mo2/(XtCO2/Mco2 + XtO2/Mo2 + XtN2/Mn2);
YtN2 = XtN2/Mn2/(XtCO2/Mco2 + XtO2/Mo2 + XtN2/Mn2);

Tt = F_T(ff,tt);
PtCO2 = (F_P(ff,tt) + Po)*YtCO2;
PtO2 = (F_P(ff,tt) + Po)*YtO2;

Rt = x1[1];
ht = (x1[1]-xc1[1])*2;

// CALCULATING FLUXES
// calculating model parameters
rho_c = 2512.8 - 0.5441*Ts;
Cc = rho_c/Mc;
Dc = (1.58E-7)*exp(-42130/(R*Ts));
Sv = Av/Ts*exp(-Ev/(R*Ts));
Sh = 100*((-4E-6)*Ts*Ts*Ts + 0.0072*Ts*Ts - 2.8306*Ts + 1330.7);

// calculating CO2 flux
term1 = e_t_MC*Cc*Dc*R*Ts*Sv;
term2 = H*(e_t_MC/e_t_SO*Zv*Zc*Cc*Dc*F*F - R*Ts*Sv);
Jco2[ii] = term1/term2*log(PtCO2/PsCO2);
Jo2[ii] = 0;

H = Rm - Rt;

fprintf(fp,"Xt=%f Xs=%f PtCO2=%f PsCO2=%f Jco2=%e \n",x1[0],x0[0],PtCO2,PsCO2,Jco2[ii]);
break;
}
}
end_f_loop(ff,tt)

fprintf(fp,"Xt=%f Xs=%f PtCO2=%f PsCO2=%f Jco2=%e \n",x1[0],x0[0],PtCO2,PsCO2,Jco2[ii]);

ii++;
}
end_f_loop(f,t)

fclose (fp);
}

[AlexanderZ]

why do you have break inside loop?
try to remove it

not sure loop over faces could be efficient inside other faces loop (check this thread, may be you may get corresponded face thread using this method UDF: Given coordinates, How to identify the cell? )

also, in case you are using fluent v19 and earlier reading/writing to file must be parallelized, check ansys fluetn customization manual

compile code

[Mohammadreza01]

Hello Alexander,
Thank you for your reply, and sorry for my late response. You helped me a lot.
I have modified my code with CX_Find_Cell_With_Point() function and it runs perfectly. I am using Fluent 2021 R1 with one processor. This is the modified code:

#include "udf.h"
#include "cxndsearch.h"
#include "surf.h"


// GEOMETRY CONSTANTS
#define L 0.015 // length of the domain [m]
#define N 170 // number of mesh cells along the domain
#define H 0.75E-3 // membrane thickness [m]

// BOUNDARY ID
#define IDmss 13 // boundary ID for membrane_shell_side
#define IDmts 12 // boundary ID for membrane_tube_side

// SPECIES ID
#define idCO2 0
#define idO2 1
#define idN2 2

// MODEL CONSTANTS
#define Po 101325 // operating pressure [Pa]
#define e_t_MC 0.145
#define e_t_SO 0.48
#define Zv 2 // oxygen vacancy charge
#define Zc -2 // carbonate ion charge
#define Zh 1 // electron hole charge
#define F 9.649E4 // Faraday's constant [C/mol]
#define R 8.3145 // ideal gas constant [J/mol-K]
#define Mc 0.02004 // molecular weight of the carbonate mixture [kg/mol]
#define Av 3.795E8 // [K/Ohm-m]
#define Ev 102415 // [J/mol]

// OTHER CONSTANTS
#define PI 3.14159265358979
#define Mco2 0.04400995 // molecular weight of CO2 [kg/mol]
#define Mo2 0.0319988 // molecular weight of O2 [kg/mol]
#define Mn2 0.0280134 // molecular weight of N2 [kg/mol]

// VARIABLES
static real Jco2[N]; // CO2 flux through membrane
static real Jo2[N]; // O2 flux through membrane
static real XmCO2; // mass fraction of CO2 in permeate
static real XmO2; // mass fraction of O2 in permeate
static real Rt; // tube radius [m]
static real Rm; // membrane radius [m]
static real hs; // boundary cell height at the shell side
static real ht; // boundary cell height at the tube side
static real Ts; // temperature at shell side
static real PsCO2; // partial pressure of CO2 at shell side
static real PsO2; // partial pressure of O2 at shell side
static real Tt; // temperature at tube side
static real PtCO2; // partial pressure of CO2 at tube side
static real PtO2; // partial pressure of O2 at tube side




DEFINE_ADJUST(flux_calculator, domain)
{
Thread *ts,*tt,*t0,*t1;
face_t f,ff;
cell_t c,c0,c1;
domain = Get_Domain(1);
ts = Lookup_Thread(domain,IDmss);
tt = Lookup_Thread(domain,IDmts);

real x0[ND_ND],xc0[ND_ND],x1[ND_ND],xc1[ND_ND],P[ND_ND];
real XsCO2,XsO2,XsN2,YsCO2,YsO2,YsN2;
real XtCO2,XtO2,XtN2,YtCO2,YtO2,YtN2;
real Cc; // carbonate ion concentration [mol/m^3]
real rho_c; // density of the carbonate mixture [kg/m^3]
real Dc; // carbonate ion diffusivity [m^2/s]
real Sv; // oxygen vacancy conductivity [1/Ohm-m]
real Sh; // electron hole conductivity [1/Ohm-m]
real term1,term2;


CX_Cell_Id *cx_cell; /* Define the current cell */
ND_Search *domain_table = NULL;
domain_table = CX_Start_ND_Point_Search(domain_table,TRUE,-1);


// GETTING DATA FROM SHELL SIDE
FILE *fp;
fp = fopen ("Output01.txt","a");

int ii = 0;
begin_f_loop(f,ts)
{
c0 = F_C0(f,ts);
t0 = THREAD_T0(ts);
C_CENTROID(xc0,c0,t0);
F_CENTROID(x0,f,ts);

// mass fractions at the shell side
XsCO2 = C_YI(c0,t0,idCO2);
XsO2 = C_YI(c0,t0,idO2);
XsN2 = C_YI(c0,t0,idN2);

// calculating mole fractions at the shell side
YsCO2 = XsCO2/Mco2/(XsCO2/Mco2 + XsO2/Mo2 + XsN2/Mn2);
YsO2 = XsO2/Mo2/(XsCO2/Mco2 + XsO2/Mo2 + XsN2/Mn2);
YsN2 = XsN2/Mn2/(XsCO2/Mco2 + XsO2/Mo2 + XsN2/Mn2);

Ts = F_T(f,ts);
PsCO2 = (F_P(f,ts) + Po)*YsCO2;
PsO2 = (F_P(f,ts) + Po)*YsO2;

Rm = x0[1];
hs = (xc0[1]-x0[1])*2;



P[0] = x0[0]; // position of cell on the tube side
P[1] = x0[1] - H - hs/2;
P[2] = x0[2];
cx_cell = CX_Find_Cell_With_Point(domain_table,P,0.0);

if(cx_cell)
{
c1 = RP_CELL(cx_cell); /* Get the marked cell */
t1 = RP_THREAD(cx_cell); /* Get the cell thread of this marked cell */

C_CENTROID(xc1,c1,t1);

// mass fractions at the tube side
XtCO2 = C_YI(c1,t1,idCO2);
XtO2 = C_YI(c1,t1,idO2);
XtN2 = C_YI(c1,t1,idN2);

// calculating mole fractions at the tube side
YtCO2 = XtCO2/Mco2/(XtCO2/Mco2 + XtO2/Mo2 + XtN2/Mn2);
YtO2 = XtO2/Mo2/(XtCO2/Mco2 + XtO2/Mo2 + XtN2/Mn2);
YtN2 = XtN2/Mn2/(XtCO2/Mco2 + XtO2/Mo2 + XtN2/Mn2);

Tt = C_T(c1,t1);
PtCO2 = (C_P(c1,t1) + Po)*YtCO2;
PtO2 = (C_P(c1,t1) + Po)*YtO2;

Rt = x1[1];
ht = (x1[1]-xc1[1])*2;

// CALCULATING FLUXES
// calculating model parameters
rho_c = 2512.8 - 0.5441*Ts;
Cc = rho_c/Mc;
Dc = (1.58E-7)*exp(-42130/(R*Ts));
Sv = Av/Ts*exp(-Ev/(R*Ts));
Sh = 100*((-4E-6)*Ts*Ts*Ts + 0.0072*Ts*Ts - 2.8306*Ts + 1330.7);

// calculating CO2 flux
term1 = e_t_MC*Cc*Dc*R*Ts*Sv;
term2 = H*(e_t_MC/e_t_SO*Zv*Zc*Cc*Dc*F*F - R*Ts*Sv);
Jco2[ii] = term1/term2*log(PtCO2/PsCO2);
Jo2[ii] = 0;
}



fprintf(fp,"Xt=%f Xs=%f PtCO2=%f PsCO2=%f Jco2=%e \n",xc1[0],x0[0],PtCO2,PsCO2,Jco2[ii]);

ii++;
}
end_f_loop(f,t)

fclose (fp);

domain_table = CX_End_ND_Point_Search(domain_table);
}