[c++]代码库
#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;
}
会python真的可以为所欲为0(回)
248天前
这里还有人吗1(回)
434天前
每天面对着电脑屏幕，敲打键盘。我所面对的并不只是代码，而是一种生活方式。0(回)
658天前
鸽子
0(回)
742天前