#include <stdio.h> |
#include <math.h> |
using namespace std; |
#include<iostream> |
#include<fstream> |
#include<stdlib.h> |
#include<ctype.h> |
#include<math.h> |
#include <cstdio> |
#include <cmath> |
#include <cstdlib> |
/*=========================================================== |
定义各参数 |
===========================================================*/ |
#define T0 353.15 |
double I ; |
#define P_in_a 1.0 //阳极进气压力 atm |
#define P_in_c 1.0 //阴极进气压力 atm |
#define ST_a 3.0 //阳极化学计量比 |
#define ST_c 3.0 //阴极化学计量比 |
#define RH_a_in 1.0 //阳极进气相对湿度 |
#define RH_c_in 1.0 //阴极相对湿度 |
#define A_in 1e-6 //流道进气截面积 |
#define A_act 25e-4 //反应活化面积 |
#define A_gdl 25e-4 //GDL层面积5cm*5cm |
#define A_c 12.5e-4 //沟道面积沟脊比为1 |
#define Sh 4.86 //舍伍德数 |
#define lm_a 0.2 //阳极CL层电解质体积分数 |
#define lm_c 0.2 //阴极CL层电解质体积分数 |
#define K_MEM 2.0e-18 //交换膜的液态水渗透率 |
#define rho_m 1980 //交换膜干态质量kg.m^-3 |
#define EW 1.1 //交换膜当量质量kg.mol^-1 |
#define alpha 0.5 //电化学传输系数 |
#define C_H2_ref 5 //氢气参考浓度 |
#define C_O2_ref 1 //氧气参考浓度 |
#define P0 101325.0 //标准大气压 |
#define Farad 96487.0 //法拉第常数 |
#define R_gas 8.314 //气体常数 |
#define Rho_l 1000.0 //液态水密度 |
#define MW_w 18e-3 //水分子摩尔质量 |
#define M_PI 3.14159265 //圆周率 |
#define Sigma_s 1500.0 //多孔层电子电导率S.m^-1 |
#define Sigma_bp 20000.0 //双极板电子电导率 |
#define Mthick 125.4e-6 //质子交换膜厚度 |
#define Gthick 300.0e-6 //GDL厚度(气体扩散层) |
#define mthick 40.0e-6 //MPL厚度(微孔层) |
#define Cthick 10.0e-6 //CL厚度(催化层) |
#define act_a 0.8 //阳极水蒸气火活度 |
#define act_c 1.0 //阴极水蒸气活度 |
#define J_l 2.0 //极限电流密度 |
#define T_c 350e-6//阴极厚度 |
#define T_a 350e-6//阳极厚度 |
double thick[4] = { 3.0e-4,0.4e-4,0.1e-4,0.5e-4 }; //各层厚度,依次为:GDL,MPL,CL和膜 |
double por[3] = { 0.6,0.4,0.3 }; //各层孔隙率,以此为GDL,MPL,CL |
double theta[3] = { 110,110,100 }; //各层接触角 |
double K0[3] = { 1.0e-12,2.5e-13,1.0e-13 }; //各层渗透率 |
//阴阳极各层液态体积分数迭代变量,依次为GDL,MPL,CL,CL-PEM |
double s_a[4] = { 0.1,0.1,0.1,0.1 }; |
double s_c[4] = { 0.5,0.5,0.5,0.5 }; |
double dT = T0 - 273.15; |
double exponent = -2.1794 + 0.02953 * dT - 9.1837e-5 * dT * dT + 0.4454e-7 * dT * dT * dT; |
double p_sat = pow (10.0, exponent) * P0; //饱和蒸汽压 |
double c_sat = p_sat / R_gas / T0; //饱和水蒸气浓度 |
//阴阳极各层水蒸气浓度迭代变量,依次为GDL,MPL,CL,CL-PEM |
double c_vap_a[4] = { c_sat, c_sat, c_sat, c_sat }; |
double c_vap_c[4] = { c_sat, c_sat, c_sat, c_sat }; |
double afax = ( double )2.25; |
double dac = ( double )2.00; |
double Get_p_sat( double T) { |
double p_sat; |
double dT = T - 273.15; |
double exponent = -2.1794 + 0.02953 * dT - 9.1837e-5 * dT * dT + 1.4454e-7 * dT * dT * dT; |
p_sat = pow (10.0, exponent) * P0; |
return p_sat; |
} |
double Get_v_thermo() { |
double v_thermo; |
v_thermo = 1.229 - 0.846e-3 * (T0 - 298.15) + R_gas * T0 / 2 / Farad * ( log (1 - Get_p_sat(T0) / P0) + 0.5 * log (0.21 * (1 - Get_p_sat(T0) / P0))); |
return v_thermo; |
} |
double water_content( double act) |
{ |
double lam = 0.0; |
if (act >= 0.0 && act <= 1.0) |
lam = 0.0043 + 17.81 * act - 39.85 * act * act + 36.0 * act * act * act; |
if (act >= 1.0 && act <= 3.0) |
lam = 14 + 1.4 * (act - 1); |
if (act > 3) { |
cout << "water content:act is big" << endl; |
lam = 16.8; |
} |
if (act < 0) { |
cout << "water content:act is small" << endl; |
lam = 0.043; |
} |
//cout << lam<<endl ; |
return lam; |
|
} |
//计算膜态水扩散系数 |
double Dm_a() |
{ |
double Dm = 0.0; |
//RH_c(2) + 2 * s_c[2]; |
double lam_a = water_content(0.8); |
double lam_c = water_content(1.0); |
double lam = (lam_a + lam_c) / 2; |
if (lam < 0.0) { |
cout << "Dm:lam is small" << endl; |
Dm = 2.693e-10; |
} |
else if (lam <= 2.0 && lam >= 0.0) |
{ |
Dm = 2.693e-10; |
} |
else if (lam > 2.0 && lam <= 3.0) |
{ |
Dm = 1.0e-10 * exp (2416 * (1.0 / 303.15 - 1.0 / T0)) * (0.87 * (3 - lam) + 2.95 * (lam - 2)); |
} |
else if (lam > 3.0 && lam <= 4.0) |
{ |
Dm = 1.0e-10 * exp (2416 * (1.0 / 303.15 - 1.0 / T0)) * (2.95 * (4 - lam) + 1.642 * (lam - 3)); |
} |
else if (lam > 4.0) |
//Dm = 1.0e-10 * exp(2416 * (1.0 / 303.15 - 1.0 / T0)) * (2.563 - 0.33 * lam + 0.0264 * lam * lam - 0.000671 * lam * lam * lam);; |
{ |
Dm = 1.0e-6 * exp (2416 * (1.0 / 303.15 - 1.0 / T0)) * (2.563 - 0.33 * lam + 0.0264 * lam * lam - 0.000671 * lam * lam * lam); |
} |
//cout << lam_a << endl << lam_c << endl; |
//cout << Dm << endl; |
return Dm; |
} |
double ndrag = 2.5; |
double sigema( double z ) //膜中水分布计算公式 |
{ |
//double afax = 1.12; |
//double dac = 2.30; |
double h_1 = ( double )(1268.00 / 303.00); |
double h_2 =( double )(1268.00 / 353.00); |
double h = h_1 - h_2; |
double ndrag = 2.5; |
//double *p ;//输入电流值 |
// p = &j; |
//double g = j * z * 1.0 * ndrag / 22 / Farad / 0.00197 / 3.81 * 10e-6; //Dm_a(); |
// double j = 0.7; |
double g = z* I * 1.0 * ndrag / (22 * Farad * 0.00197 * Dm_a()); |
double lamz = (0.005193*(4.4 * afax + dac * exp (g))-0.00326)* exp (h); |
// cout << "g="<< g << endl; |
//cout << "h=" << h<< endl; |
return lamz; |
|
} |
double Rm( double z) |
{ |
double Rm_a = 1 / sigema(z); |
return Rm_a; |
} |
double f( double z) |
{ |
return Rm(z); |
} |
double adaptive_simpson( double a, double b, double d) |
{ |
double c = a + (b - a) / 2; |
double l = (f(a) + 4 * f(c) + f(b)) * (b - a) / 6; |
double r = (f(a) + 4 * f((a + c) / 2) + 2 * f(c) + 4 * f((c + b) / 2) + f(b)) * (b - a) / 12; |
if ( fabs (l - r) <= 15 * d) |
return r + (r - l) / 15; |
return adaptive_simpson(a, c, d) + adaptive_simpson(c, b, d); |
} |
double Rm_jf() //自适应辛普森公式算积分 |
{ |
double a = 0; //积分下限 |
double b = 0.0125; //积分上限 |
double d = 0.00000001; //精度 |
double rm_jf = adaptive_simpson(a, b, d); |
//cout <<"rm_jf=" << rm_jf << endl; |
return rm_jf; |
} |
double Get_ohmic( ) |
{ |
|
double ohm = I * Rm_jf(); |
//cout << "ohm=" << ohm << endl; |
return ohm; |
|
} |
/*double s_j() |
{ |
double j = I; |
return j; |
}*/ |
//void anotherfunction(); |
double C_cha( ) //浓度损耗 |
{ |
double cha = R_gas * T0 * (1.0 + 1 / 0.5) / (2* Farad); |
return cha; |
} |
double Get_conc() |
{ |
double conc = C_cha() * ( log (J_l)- log (J_l-I)); |
// cout << "conc=" << conc << endl; |
return conc; |
} |
/*阴极活化电势,Ω.m2*/ |
double Get_eta_c() |
{ |
double j0_a_ref = 0.001; //交换电流密度 |
double P_c = 3; //阴极压强 |
double x_o2 = 0.19; //氧气摩尔分数 |
double Dl_o = 0; |
double lam = water_content(1.0); |
if (lam < 0.0) { |
cout << "Dm:lam is small" << endl; |
Dl_o = 2.693e-10; |
} |
if (lam <= 2.0 && lam >= 0.0) |
{ |
Dl_o = 2.693e-10; |
} |
else if (lam > 2.0 && lam <= 3.0) |
{ |
Dl_o = 1.0e-10 * exp (2416 * (1.0 / 303.15 - 1.0 / T0)) * (0.87 * (3 - lam) + 2.95 * (lam - 2)); |
} |
else if (lam > 3.0 && lam <= 4.0) |
{ |
Dl_o = 1.0e-10 * exp (2416 * (1.0 / 303.15 - 1.0 / T0)) * (2.95 * (4 - lam) + 1.642 * (lam - 3)); |
} |
else if (lam > 4.0) |
{ |
Dl_o = 1.0e-6 * exp (2416 * (1.0 / 303.15 - 1.0 / T0)) * (2.563 - 0.33 * lam + 0.0264 * lam * lam - 0.000671 * lam * lam * lam); |
} |
// cout << "Dm=" << Dm << endl; |
|
double down = j0_a_ref*P_c*(x_o2- T_c*I* R_gas * T0/(4* Farad * P_c* Dl_o)); |
double eta_c = R_gas * T0 / alpha / 4 / Farad * ( log (I)- log (down)); |
return eta_c; |
} |
/*阳极活化电势,Ω.m2*/ |
double Get_eta_a() |
{ |
double j0_a_ref = 0.001; //交换电流密度 |
double P_a = 3; //阳极压强 |
double x_h2 = 0.9; //氧气摩尔分数 |
double Dl_o = 0; |
double lam = water_content(0.8); |
if (lam < 0.0) { |
cout << "Dm:lam is small" << endl; |
Dl_o = 2.693e-10; |
} |
if (lam <= 2.0 && lam >= 0.0) |
{ |
Dl_o = 2.693e-10; |
} |
else if (lam > 2.0 && lam <= 3.0) |
{ |
Dl_o = 1.0e-10 * exp (2416 * (1.0 / 303.15 - 1.0 / T0)) * (0.87 * (3 - lam) + 2.95 * (lam - 2)); |
} |
else if (lam > 3.0 && lam <= 4.0) |
{ |
Dl_o = 1.0e-10 * exp (2416 * (1.0 / 303.15 - 1.0 / T0)) * (2.95 * (4 - lam) + 1.642 * (lam - 3)); |
} |
else if (lam > 4.0) |
{ |
Dl_o = 1.0e-6 * exp (2416 * (1.0 / 303.15 - 1.0 / T0)) * (2.563 - 0.33 * lam + 0.0264 * lam * lam - 0.000671 * lam * lam * lam); |
} |
// cout << "Dm=" << Dl_o << endl; |
double down = j0_a_ref * P_a * (x_h2 - T_c * I * R_gas * T0 / (4 * Farad * P_a * Dl_o)); //需要调整 |
double eta_a = R_gas * T0 / alpha / 2 / Farad * ( log (I) - log (down)); |
return eta_a; |
} |
int main() |
{ |
|
cout << "请输入电流值,以-1结束\n" ; |
double CURRENT_DENSITY[30] = { 0 }; |
for ( int i = 0; i < 30; i++) { |
cin >> CURRENT_DENSITY[i]; |
if (CURRENT_DENSITY[i] == -1) |
break ; |
} |
int j = 0; |
while (j < 30 && CURRENT_DENSITY[j] != -1) { |
I = CURRENT_DENSITY[j]; |
//bool che = 0; |
//double I = 0.7; |
double V = Get_v_thermo() - Get_ohmic( )- Get_conc()- Get_eta_c()- Get_eta_a(); |
|
j++; |
cout << V<<endl<< Get_v_thermo()<<endl<< Get_ohmic()<<endl << Get_conc() << endl; |
cout << Get_eta_a() << endl<< Get_eta_c(); |
} |
system ( "pause" ); |
return 0; |
} |