一维仿真

[c++]代码库

#include <iostream>
#include <cmath>
#include <vector>
using namespace std;

// 定义一些常数和参数
const double F = 96485.0; // 法拉第常数 (C/mol)
const double R = 8.314; // 通用气体常数 (J/mol/K)
const double T = 353.0; // 工作温度 (K)
const double i0 = 0.1; // 交换电流密度 (A/cm2)
const double alpha = 0.5; // 电荷转移系数
const double b = 0.05; // 质量传递系数 (cm/s)
const double Lm = 25e-4; // 膜厚度 (cm)
const double Lc = 10e-4; // 催化层厚度 (cm)
const double Lg = 200e-4; // 气体扩散层厚度 (cm)
const double sigma_m = 0.1; // 膜电导率 (S/cm)

const double sigma_c = 100.0; // 催化层电导率 (S/cm)
const double sigma_g = 10.0; // 气体扩散层电导率 (S/cm)
const double cH2 = 0.2; // 氢气浓度 (mol/cm3)
const double cO2 = 0.1; // 氧气浓度 (mol/cm3)

// 定义边界条件
const double phi_a = 0.0; // 阳极电势 (V)
const double phi_c = 0.6; // 阴极电势 (V)

// 定义网格大小和节点数
const int Nm = 5; // 膜中的节点数
const int Nc = 5; // 催化层中的节点数
const int Ng = 20; // 气体扩散层中的节点数
const int N = Nm + 2 * Nc + 2 * Ng + 1; // 总节点数
const double dxm = Lm / Nm; // 膜中的网格大小
const double dxc = Lc / Nc; // 催化层中的网格大小
const double dxg = Lg / Ng; // 气体扩散层中的网格大小

// 定义一个函数来计算某个节点处的过电势
double eta(int i, double phi, double c) {
    if (i <= Ng + Nc) { // 阳极侧
        return phi - log(cH2 / c) / (2 * alpha * F / R / T) - i0 / sigma_c / dxc;
    }
    else if (i >= Ng + Nc + Nm + Nc) { // 阴极侧
        return phi - log(cO2 / c) / (4 * alpha * F / R / T) - i0 / sigma_c / dxc;
    }
    else { // 膜或界面
        return 0.0;
    }
}

// 定义一个函数来计算某个节点处的电流密度
double current(int i, double phi, double c) {
    return i0 * (exp(2 * alpha * F * eta(i, phi, c) / R / T) - exp(-2 * (1 - alpha) * F * eta(i, phi, c) / R / T));
}

// 定义一个函数来计算某个节点处的浓度
double concentration(int i, double phi, double c) {
    if (i <= Ng + Nc) { // 阳极侧
        return cH2 * exp(-b * current(i, phi, c));
    }
    else if (i >= Ng + Nc + Nm + Nc) { // 阴极侧
        return cO2 * exp(-b * current(i, phi, c));

    }
    else { // 膜或界面
        return c;
    }
}

// 定义一个函数来计算某个节点处的电导率
double conductivity(int i) {
    if (i <= Ng) { // 阳极气体扩散层
        return sigma_g;
    }
    else if (i <= Ng + Nc) { // 阳极催化层
        return sigma_c;
    }
    else if (i <= Ng + Nc + Nm) { // 膜
        return sigma_m;
    }
    else if (i <= Ng + Nc + Nm + Nc) { // 阴极催化层
        return sigma_c;
    }
    else { // 阴极气体扩散层
        return sigma_g;
    }
}

// 定义一个函数来计算某个节点处的网格大小
double mesh_size(int i) {
    if (i <= Ng) { // 阳极气体扩散层
        return dxg;
    }
    else if (i <= Ng + Nc) { // 阳极催化层
        return dxc;
    }
    else if (i <= Ng + Nc + Nm) { // 膜
        return dxm;
    }
    else if (i <= Ng + Nc + Nm + Nc) { // 阴极催化层
        return dxc;
    }
    else { // 阴极气体扩散层
        return dxg;
    }
}

// 定义一个向量来存储每个节点处的电势
vector<double> phi(N);

// 定义一个向量来存储每个节点处的浓度
vector<double> c(N);

// 定义一个向量来存储每个节点处的电流密度
vector<double> i(N);

// 用线性插值法初始化电势和浓度向量
for (int j = 0; j < N; j++) {
    phi[j] = phi_a + (phi_c - phi_a) * j / (N - 1);

    c[j] = 0.5 * (cH2 + cO2);
}

// 定义一个收敛性判断的容差
const double tol = 1e-6;

// 定义一个变量来存储最大误差
double err;

// 定义一个变量来计数迭代次数
int iter = 0;

// 开始一个循环，直到收敛或达到最大迭代次数
do {
    // 用最新的电势和浓度值更新电流密度向量
    for (int k = 0; k < N; k++) {
        i[k] = current(k, phi[k], c[k]);
    }

    // 用高斯-赛德尔法求解离散化的泊松方程，更新电势向量
    for (int l = 1; l < N - 1; l++) {
        phi[l] = (conductivity(l - 1) * phi[l - 1] + conductivity(l + 1) * phi[l + 1] - i[l] * mesh_size(l)) / (conductivity(l - 1) + conductivity(l + 1));
    }

    // 用高斯-赛德尔法求解离散化的质量守恒方程，更新浓度向量
    for (int m = 1; m < N - 1; m++) {
        c[m] = (c[m - 1] + c[m + 1]) / 2 - b * i[m] * mesh_size(m) / 2;
    }

    // 计算旧值和新值之间的最大误差
    err = 0.0;
    for (int n = 0; n < N; n++) {
        err = max(err, abs(phi[n] - phi_old[n]));
        err = max(err, abs(c[n] - c_old[n]));
    }

    // 增加迭代次数
    iter++;

    // 打印迭代次数和误差
    cout << "Iteration " << iter << ", error = " << err << endl;

} while (err > tol && iter < 1000); // 循环结束

// 检查是否收敛
if (err <= tol) {
    cout << "Convergence achieved!" << endl;
}
else {
    cout << "Convergence failed!" << endl;
}

// 打印结果
cout << "Node\tPotential\tConcentration\tCurrent" << endl;
for (int p = 0; p < N; p++) {
    cout << p << "\t" << phi[p] << "\t" << c[p] << "\t" << i[p] << endl;
}

// 定义一个主函数来运行仿真程序
int main() {
    // 调用仿真函数
    simulate();

    // 返回0表示成功执行
    return 0;
}