
[Sponsors] 
October 15, 2015, 06:02 
UDS equation for liquid water saturation in porous media

#1 
Member
Join Date: Jul 2014
Location: Portugal
Posts: 34
Rep Power: 5 
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! 

October 15, 2015, 10:16 

#2 
Senior Member
Bruno Machado
Join Date: May 2014
Posts: 271
Rep Power: 6 
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. 

October 15, 2015, 13:01 

#3 
Member
Join Date: Jul 2014
Location: Portugal
Posts: 34
Rep Power: 5 
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. 

October 15, 2015, 14:13 

#4 
Senior Member
Bruno Machado
Join Date: May 2014
Posts: 271
Rep Power: 6 
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. 

October 16, 2015, 05:29 

#5 
Member
Join Date: Jul 2014
Location: Portugal
Posts: 34
Rep Power: 5 
Bruno, I just sent you a email with a response.
If anyone else has any other suggestions, they will be greatly appreciated! 

October 20, 2015, 10:33 

#6 
Senior Member
Bruno
Join Date: Mar 2009
Location: Brazil
Posts: 279
Rep Power: 14 
Try posting the equations you're solving, boundary conditions, etc. Maybe with more info someone could help finding the problem.


October 20, 2015, 11:57 

#7 
Member
Join Date: Jul 2014
Location: Portugal
Posts: 34
Rep Power: 5 
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 nonporous 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 nonporous 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*T9.1387/100000*pow(T,2)+1.4451/10000000*pow(T,3)); dens_liquid=0.0036*pow(T,2)0.0695*T+1000.5; /*kg m3*/ 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 = (1s)*(P_h2oPsat)/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_h2oPsat)/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_h2oPsat)/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.55e13; /*m2*/ surf_tension = 0.0002*T+0.0761; /*N m2*/ visc_liquid= (0.0026*pow(T,3)+0.5874*pow(T,2)47.598*T+1763.4)*1e6; /*Pa.s*/ dens_liquid=0.0036*pow(T,2)0.0695*T+1000.5; /*kg m3*/ K_rel=s*s*s; if (cos(GDL_c_angle*2*pi/360)>0.0) ss=1s; 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=1e20; 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=1e20; } C_UDMI(c0,t0,73)=flux; return flux; } 

January 16, 2017, 01:52 

#8  
New Member
Mahbod MoeinJahromi
Join Date: Jan 2017
Posts: 2
Rep Power: 0 
Quote:
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. 

January 16, 2017, 05:55 

#9 
Member
Join Date: Jul 2014
Location: Portugal
Posts: 34
Rep Power: 5 
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 

January 16, 2017, 07:45 

#10  
Senior Member
Bruno Machado
Join Date: May 2014
Posts: 271
Rep Power: 6 
Quote:
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 Code:
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*T9.1387/100000*pow(T,2)+1.4451/10000000*pow(T,3)); dens_liquid=0.0036*pow(T,2)0.0695*T+1000.5; /*kg m3*/ 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 = (1s)*(P_h2oPsat)/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_h2oPsat)/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_h2oPsat)/R/(T+273.15)*(MH2O*0.001)*kc; return source; } 

February 8, 2017, 11:27 

#11 
Member
Join Date: Jul 2014
Location: Portugal
Posts: 34
Rep Power: 5 
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 

September 19, 2017, 20:44 

#12  
New Member
Arif
Join Date: Sep 2017
Posts: 12
Rep Power: 2 
Quote:
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 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
suGWFoam: Richards equation solver for porous media flows  liu  OpenFOAM Announcements from Other Sources  0  November 11, 2012 21:09 
UDS for Porous Media  cp  FLUENT  6  June 26, 2012 13:44 
Multiphase Porous Media Flow  Convergence Issues  VT_Bromley  FLUENT  5  March 23, 2011 11:38 
debug of a UDS about porous media!  whiz  FLUENT  0  February 16, 2009 10:41 
porous media: Fluent or StarCD?  Igor  Main CFD Forum  0  December 5, 2002 16:16 