[Rui_27]

Hello,

I am simulating the cathode of a PEM fuel cell. I wrote a udf which accounts for the reaction of O2 into the catalyst layer and the consequent formation of water vapor.
I am now trying to account for the water vapor condensation by solving a UDS equation for the liquid water saturation in the porous media. I managed to write the UDFs for the UDS diffusivity, flux and source, but i am not getting promise results. The solution either diverges rapidly or, if not, the results are senseless.
If fact, I am trying to do the same as in the PEMFC module from Fluent, but without using their module.

Does anyone have experience on liquid water saturation simulation in porous media, and has some suggestions? Any comment will be appreciated.

Thanks!

[Bruno Machado]

I am building an entire fuel cell module as an UDF. I was getting problem in my liquid equation due to the initial condition. Depending on the way you set your initial condition, evaporation/condensation can happen and it leads to a negative liquid fraction.

Have a look at it and let me know if it worked for you.

By the way, how did you built your mesh? I am facing problem with mine =/

Regards.

[Rui_27]

Hi Bruno (Português?)

I also did notice that the initial value for the uds (liquid water saturation) is important, and I tried several options. The results did change but they are still far from acceptable (water distribution is senseless and saturation in some areas is greater than 1.0..)
Did follow any specific criteria for choosing the initial value?

Do you also work with PEM fuel cell type?

I built my mesh using the Design Modeler (Geometry) and the ANSYS Meshing (mesh) and so far so good.

Regards.

[Bruno Machado]

Hi Rui_27,

Im from Brazil, but I am doing my PhD in the UK.

I did not follow any specific criteria, but I considered 0.1 as the initial value for the liquid water and it worked fine (I tried less than this and I got problem and senseless values). In addition to that, I add a if condition to my code that if the liquid saturation scalar is equal to 0.00001, the source term is equal to 0 (no condensation/evaporation), so the liquid phase remains positive.

Yes, I am working with PEM fuel cell, but I am studying a Anion Exchange Membrane (AEM) fuel cell.

Could you please share how you did your mesh? I am not sure how to do regarding the parts/bodys and also the contact regions. Or if it wont cause you any harm, could you share your mesh/geometry file?

I will let my email here, so you can speak in Portuguese and once we solve our issues, we can post the solution here.

My email is: brunosouzamachado@gmail.com

All the best.

[Rui_27]

Bruno, I just sent you a email with a response.

If anyone else has any other suggestions, they will be greatly appreciated!

[brunoc]

Try posting the equations you're solving, boundary conditions, etc. Maybe with more info someone could help finding the problem.

[Rui_27]

Good ideia Brunoc.

I am trying to solve a UDS equation for the liquid water saturation in a cathode of a PEM fuel cell.

- First I have two UDFs with sources terms that account for the condensation/evaporation of water. One applied to the equation of h2o (vapor) and other (same value with opposite sign) applied in the UDS equation (liquid water).
- Then I have a UDF for the diffusivity of the UDS. In the porous media it corresponds to the capillary difusivity and in the non-porous media I set the diffusivity as 0.0.
- I also have a UDF for the UDS flux. It is equal to dens_liquid*gas_velocity in the non-porous zone and equal to 0.0 in the porous media.

The UDFs codes are below.

Equations are in the file attached (1.26 for nonporous and 1.28 for porous media).

Regards

-----------------------------

DEFINE_SOURCE(liquid_source_h2o,c,t,dS,eqn)
{
real P_h2o;
real Psat;
real dens_liquid;
real s;
real term1;
real term2;
real source;

s =MAX(0.0001, MIN(C_UDSI(c,t,0),0.999));
Psat= 100000.0*pow(10,-2.1794+0.02953*T-9.1387/100000*pow(T,2)+1.4451/10000000*pow(T,3));
dens_liquid=-0.0036*pow(T,2)-0.0695*T+1000.5; /*kg m-3*/
P_h2o= C_YI(c, t, 1)*C_R(c,t)*R*(T+273.15)/(MH2O*0.001);
C_UDMI(c,t,66) =P_h2o;

term1 = (1-s)*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001);
term2 = -1.0*s*dens_liquid;
C_UDMI(c,t,67) =term1;
C_UDMI(c,t,68) =term2;

if (term1>term2)
{
source = -1.0*kc*term1;
dS[eqn] = 1.0*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001)*kc;
}
else
{
source = -1.0*kc*term2;
dS[eqn] = 1.0*dens_liquid*kc;
}

C_UDMI(c,t,69) =source;
C_UDMI(c,t,70) =dS[eqn];
C_UDMI(c,t,71) = -1.0*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001)*kc;

return source;

}

DEFINE_SOURCE(liquid_source_s,c,t,dS,eqn)
{
real source;
source= -1.0*C_UDMI(c,t,69);
dS[eqn]=-1.0*C_UDMI(c,t,70);
C_UDMI(c,t,81) =source;
C_UDMI(c,t,82) =dS[eqn];

return source;
}

DEFINE_DIFFUSIVITY(liquid_diffusivity2, c, t, i)
{
real K_abs;
real surf_tension;
real visc_liquid;
real dens_liquid;
real dpcds=0.0;
real K_rel=0.0;
real s;
real ss=0.0;
real diff;
s = MAX(0.0001, MIN(C_UDSI(c,t,0),0.999));

K_abs = 2.55e-13; /*m2*/
surf_tension = -0.0002*T+0.0761; /*N m-2*/
visc_liquid= (-0.0026*pow(T,3)+0.5874*pow(T,2)-47.598*T+1763.4)*1e-6; /*Pa.s*/
dens_liquid=-0.0036*pow(T,2)-0.0695*T+1000.5; /*kg m-3*/
K_rel=s*s*s;

if (cos(GDL_c_angle*2*pi/360)>0.0)
ss=1-s;
else
ss=s;

dpcds = surf_tension*ABS(cos(GDL_c_angle*2*pi/360))/sqrt(K_abs/C_POR(c, t))*(1.417 + 3*1.263*ss*ss - 2*2.12*ss);

if (C_POR(c,t)==1.0)
diff=1e-20;
else
diff= dens_liquid*K_abs*K_rel*dpcds/visc_liquid;

C_UDMI(c,t,80)=diff;

return diff;
}

DEFINE_UDS_FLUX(uds_flux2,f,t,i)
{

cell_t c0, c1 = -1;
Thread *t0, *t1 = NULL;
real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0;
c0 = F_C0(f,t);
t0 = F_C0_THREAD(f,t);
F_AREA(A, f, t);


if (C_POR(c0,t0)==1.0)
{
if (BOUNDARY_FACE_THREAD_P(t))
{
real dens;

if (NNULLP(THREAD_STORAGE(t,SV_DENSITY)))
dens = (-0.0036*pow(T,2)-0.0695*T+1000.5);
else
dens = (-0.0036*pow(T,2)-0.0695*T+1000.5);
NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, dens);
flux = NV_DOT(psi_vec, A);
}
else
{
c1 = F_C1(f,t);
t1 = F_C1_THREAD(f,t);
NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,dens);
NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,dens);
flux = NV_DOT(psi_vec, A)/2.0;
}
}
else
{
flux=1e-20;
}

C_UDMI(c0,t0,73)=flux;

return flux;
}

[mahbod]

Originally Posted by Rui_27

Good ideia Brunoc.

I am trying to solve a UDS equation for the liquid water saturation in a cathode of a PEM fuel cell.

- First I have two UDFs with sources terms that account for the condensation/evaporation of water. One applied to the equation of h2o (vapor) and other (same value with opposite sign) applied in the UDS equation (liquid water).
- Then I have a UDF for the diffusivity of the UDS. In the porous media it corresponds to the capillary difusivity and in the non-porous media I set the diffusivity as 0.0.
- I also have a UDF for the UDS flux. It is equal to dens_liquid*gas_velocity in the non-porous zone and equal to 0.0 in the porous media.

The UDFs codes are below.

Equations are in the file attached (1.26 for nonporous and 1.28 for porous media).

Regards

-----------------------------

DEFINE_SOURCE(liquid_source_h2o,c,t,dS,eqn)
{
real P_h2o;
real Psat;
real dens_liquid;
real s;
real term1;
real term2;
real source;

s =MAX(0.0001, MIN(C_UDSI(c,t,0),0.999));
Psat= 100000.0*pow(10,-2.1794+0.02953*T-9.1387/100000*pow(T,2)+1.4451/10000000*pow(T,3));
dens_liquid=-0.0036*pow(T,2)-0.0695*T+1000.5; /*kg m-3*/
P_h2o= C_YI(c, t, 1)*C_R(c,t)*R*(T+273.15)/(MH2O*0.001);
C_UDMI(c,t,66) =P_h2o;

term1 = (1-s)*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001);
term2 = -1.0*s*dens_liquid;
C_UDMI(c,t,67) =term1;
C_UDMI(c,t,68) =term2;

if (term1>term2)
{
source = -1.0*kc*term1;
dS[eqn] = 1.0*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001)*kc;
}
else
{
source = -1.0*kc*term2;
dS[eqn] = 1.0*dens_liquid*kc;
}

C_UDMI(c,t,69) =source;
C_UDMI(c,t,70) =dS[eqn];
C_UDMI(c,t,71) = -1.0*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001)*kc;

return source;

}

DEFINE_SOURCE(liquid_source_s,c,t,dS,eqn)
{
real source;
source= -1.0*C_UDMI(c,t,69);
dS[eqn]=-1.0*C_UDMI(c,t,70);
C_UDMI(c,t,81) =source;
C_UDMI(c,t,82) =dS[eqn];

return source;
}

DEFINE_DIFFUSIVITY(liquid_diffusivity2, c, t, i)
{
real K_abs;
real surf_tension;
real visc_liquid;
real dens_liquid;
real dpcds=0.0;
real K_rel=0.0;
real s;
real ss=0.0;
real diff;
s = MAX(0.0001, MIN(C_UDSI(c,t,0),0.999));

K_abs = 2.55e-13; /*m2*/
surf_tension = -0.0002*T+0.0761; /*N m-2*/
visc_liquid= (-0.0026*pow(T,3)+0.5874*pow(T,2)-47.598*T+1763.4)*1e-6; /*Pa.s*/
dens_liquid=-0.0036*pow(T,2)-0.0695*T+1000.5; /*kg m-3*/
K_rel=s*s*s;

if (cos(GDL_c_angle*2*pi/360)>0.0)
ss=1-s;
else
ss=s;

dpcds = surf_tension*ABS(cos(GDL_c_angle*2*pi/360))/sqrt(K_abs/C_POR(c, t))*(1.417 + 3*1.263*ss*ss - 2*2.12*ss);

if (C_POR(c,t)==1.0)
diff=1e-20;
else
diff= dens_liquid*K_abs*K_rel*dpcds/visc_liquid;

C_UDMI(c,t,80)=diff;

return diff;
}

DEFINE_UDS_FLUX(uds_flux2,f,t,i)
{

cell_t c0, c1 = -1;
Thread *t0, *t1 = NULL;
real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0;
c0 = F_C0(f,t);
t0 = F_C0_THREAD(f,t);
F_AREA(A, f, t);


if (C_POR(c0,t0)==1.0)
{
if (BOUNDARY_FACE_THREAD_P(t))
{
real dens;

if (NNULLP(THREAD_STORAGE(t,SV_DENSITY)))
dens = (-0.0036*pow(T,2)-0.0695*T+1000.5);
else
dens = (-0.0036*pow(T,2)-0.0695*T+1000.5);
NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, dens);
flux = NV_DOT(psi_vec, A);
}
else
{
c1 = F_C1(f,t);
t1 = F_C1_THREAD(f,t);
NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,dens);
NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,dens);
flux = NV_DOT(psi_vec, A)/2.0;
}
}
else
{
flux=1e-20;
}

C_UDMI(c0,t0,73)=flux;

return flux;
}

Hi, Rui_27
Actually, I faced the same problem as yours about the UDS equation for liquid water saturation in porous media .
I followed this topic in order to find my answer, but after you udf/uds code there is nothing. Did you solve your problem? How? Please give me a hint to solve my problem thanks.

[Rui_27]

I Mahbod,

Yes. I "solved" it. After a lot of trial and error I found that setting the minimal value of s equal to 0.02 in DEFINE_DIFFUSIVITY(liquid_diffusivity2, c, t, i), instead of 0.0001, solves my problem.

Liquid diffusivity changes greatly with small variations of s, specially at low values of s, and this causes difficulties in converging the UDS equations. Preventing s from getting to small, in my case, helped.

The other way to solve this is to change the solution controls apropriately (under relaxation factors, cycle types, ect.), but I did not explore this strategy in detail. See more on this in:

Hao Wu. Mathematical modeling of transient transport phenomena in PEM fuel cells. Canada: University of Waterloo; 2009.

Regards,
Rui

[Bruno Machado]

Originally Posted by Rui_27

I Mahbod,

Yes. I "solved" it. After a lot of trial and error I found that setting the minimal value of s equal to 0.02 in DEFINE_DIFFUSIVITY(liquid_diffusivity2, c, t, i), instead of 0.0001, solves my problem.

Liquid diffusivity changes greatly with small variations of s, specially at low values of s, and this causes difficulties in converging the UDS equations. Preventing s from getting to small, in my case, helped.

The other way to solve this is to change the solution controls apropriately (under relaxation factors, cycle types, ect.), but I did not explore this strategy in detail. See more on this in:

Hao Wu. Mathematical modeling of transient transport phenomena in PEM fuel cells. Canada: University of Waterloo; 2009.

Regards,
Rui

Hi Rui,

I hope to find you well. I am glad that you make your code work. Well done.

Just a question, are you still using this for condensation AND evaporation? I am asking because in the manual you posted, it only defines the condensation, and I was wondering where did you see this approach for evaporation.

Regards,
Bruno

 DEFINE_SOURCE(liquid_source_h2o,c,t,dS,eqn)
{
real P_h2o;
real Psat;
real dens_liquid;
real s;
real term1;
real term2;
real source;

s =MAX(0.0001, MIN(C_UDSI(c,t,0),0.999));
Psat= 100000.0*pow(10,-2.1794+0.02953*T-9.1387/100000*pow(T,2)+1.4451/10000000*pow(T,3)); 
dens_liquid=-0.0036*pow(T,2)-0.0695*T+1000.5; /*kg m-3*/ 
P_h2o= C_YI(c, t, 1)*C_R(c,t)*R*(T+273.15)/(MH2O*0.001);
C_UDMI(c,t,66) =P_h2o;

term1 = (1-s)*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001);
term2 = -1.0*s*dens_liquid;
C_UDMI(c,t,67) =term1;
C_UDMI(c,t,68) =term2;

if (term1>term2)
{
source = -1.0*kc*term1;
dS[eqn] = 1.0*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001)*kc;
}
else
{
source = -1.0*kc*term2;
dS[eqn] = 1.0*dens_liquid*kc;
} 

C_UDMI(c,t,69) =source;
C_UDMI(c,t,70) =dS[eqn];
C_UDMI(c,t,71) = -1.0*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001)*kc;

return source;

}

[Rui_27]

Hi Bruno!

Sorry for the late response.

Yes. When the Source in MACRO liquid_source_h2o is negative, there is condensation. When it is positive, there is evaporation.

I called the MACRO liquid_source_h2o, but in fact in this macro computes the source term for water vapor. The source term in MACRO liquid_source_s is set as the mass source term for liquid water.

In my model i did not consider heat transfer, so condensation/evaporation is just due to local concentration variations.

Regards,
Rui

[engrarif]

Originally Posted by Rui_27

Good ideia Brunoc.

I am trying to solve a UDS equation for the liquid water saturation in a cathode of a PEM fuel cell.

- First I have two UDFs with sources terms that account for the condensation/evaporation of water. One applied to the equation of h2o (vapor) and other (same value with opposite sign) applied in the UDS equation (liquid water).
- Then I have a UDF for the diffusivity of the UDS. In the porous media it corresponds to the capillary difusivity and in the non-porous media I set the diffusivity as 0.0.
- I also have a UDF for the UDS flux. It is equal to dens_liquid*gas_velocity in the non-porous zone and equal to 0.0 in the porous media.

The UDFs codes are below.

Equations are in the file attached (1.26 for nonporous and 1.28 for porous media).

Regards

-----------------------------

DEFINE_SOURCE(liquid_source_h2o,c,t,dS,eqn)
{
real P_h2o;
real Psat;
real dens_liquid;
real s;
real term1;
real term2;
real source;

s =MAX(0.0001, MIN(C_UDSI(c,t,0),0.999));
Psat= 100000.0*pow(10,-2.1794+0.02953*T-9.1387/100000*pow(T,2)+1.4451/10000000*pow(T,3));
dens_liquid=-0.0036*pow(T,2)-0.0695*T+1000.5; /*kg m-3*/
P_h2o= C_YI(c, t, 1)*C_R(c,t)*R*(T+273.15)/(MH2O*0.001);
C_UDMI(c,t,66) =P_h2o;

term1 = (1-s)*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001);
term2 = -1.0*s*dens_liquid;
C_UDMI(c,t,67) =term1;
C_UDMI(c,t,68) =term2;

if (term1>term2)
{
source = -1.0*kc*term1;
dS[eqn] = 1.0*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001)*kc;
}
else
{
source = -1.0*kc*term2;
dS[eqn] = 1.0*dens_liquid*kc;
}

C_UDMI(c,t,69) =source;
C_UDMI(c,t,70) =dS[eqn];
C_UDMI(c,t,71) = -1.0*(P_h2o-Psat)/R/(T+273.15)*(MH2O*0.001)*kc;

return source;

}

DEFINE_SOURCE(liquid_source_s,c,t,dS,eqn)
{
real source;
source= -1.0*C_UDMI(c,t,69);
dS[eqn]=-1.0*C_UDMI(c,t,70);
C_UDMI(c,t,81) =source;
C_UDMI(c,t,82) =dS[eqn];

return source;
}

DEFINE_DIFFUSIVITY(liquid_diffusivity2, c, t, i)
{
real K_abs;
real surf_tension;
real visc_liquid;
real dens_liquid;
real dpcds=0.0;
real K_rel=0.0;
real s;
real ss=0.0;
real diff;
s = MAX(0.0001, MIN(C_UDSI(c,t,0),0.999));

K_abs = 2.55e-13; /*m2*/
surf_tension = -0.0002*T+0.0761; /*N m-2*/
visc_liquid= (-0.0026*pow(T,3)+0.5874*pow(T,2)-47.598*T+1763.4)*1e-6; /*Pa.s*/
dens_liquid=-0.0036*pow(T,2)-0.0695*T+1000.5; /*kg m-3*/
K_rel=s*s*s;

if (cos(GDL_c_angle*2*pi/360)>0.0)
ss=1-s;
else
ss=s;

dpcds = surf_tension*ABS(cos(GDL_c_angle*2*pi/360))/sqrt(K_abs/C_POR(c, t))*(1.417 + 3*1.263*ss*ss - 2*2.12*ss);

if (C_POR(c,t)==1.0)
diff=1e-20;
else
diff= dens_liquid*K_abs*K_rel*dpcds/visc_liquid;

C_UDMI(c,t,80)=diff;

return diff;
}

DEFINE_UDS_FLUX(uds_flux2,f,t,i)
{

cell_t c0, c1 = -1;
Thread *t0, *t1 = NULL;
real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0;
c0 = F_C0(f,t);
t0 = F_C0_THREAD(f,t);
F_AREA(A, f, t);


if (C_POR(c0,t0)==1.0)
{
if (BOUNDARY_FACE_THREAD_P(t))
{
real dens;

if (NNULLP(THREAD_STORAGE(t,SV_DENSITY)))
dens = (-0.0036*pow(T,2)-0.0695*T+1000.5);
else
dens = (-0.0036*pow(T,2)-0.0695*T+1000.5);
NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, dens);
flux = NV_DOT(psi_vec, A);
}
else
{
c1 = F_C1(f,t);
t1 = F_C1_THREAD(f,t);
NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,dens);
NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,dens);
flux = NV_DOT(psi_vec, A)/2.0;
}
}
else
{
flux=1e-20;
}

C_UDMI(c0,t0,73)=flux;

return flux;
}

Hi, RUI,
Great work, i really appreciate it. I have a question that how can we define the value of kc (condensation rate constant) ??? i just need to change this value from 100/s to 500/s , as it is hardwired by fluent to a value of 100/ s
Thanks in advance