一维仿真-1

[c++]代码库

#include <iostream>
#include <cmath>
#include <cstdlib>
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 dx(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;
    }
}

// 定义一个函数来使用追赶法解线性方程组 Ax = b，其中 A 是三对角矩阵，x 和 b 是向量，n 是维数
void thomas(double* a, double* b, double* c, double* d, double* x, int n) {
    // a, b, c 是三对角矩阵 A 的对角线元素，d 是右端向量 b，x 是解向量 x
    // 使用一维数组存储三对角矩阵和向量，下标从 1 开始
    double* p = new double[n]; // 用于存储中间变量 p
    double* q = new double[n]; // 用于存储中间变量 q

    // 前向消元过程，得到 p 和 q
    p[1] = b[1];
    q[1] = d[1] / p[1];
    for (int i = 2; i < n; i++) {
        p[i] = b[i] - a[i] * c[i - 1] / p[i - 1];
        q[i] = (d[i] - a[i] * q[i - 1]) / p[i];
    }
    
    // 后向回代过程，得到 x
    x[n - 1] = q[n - 1];
    for (int i = n - 2; i >= 1; i--) {
        x[i] = q[i] - c[i] * x[i + 1] / p[i];
    }

    // 释放动态分配的内存空间
    delete[] p;
    delete[] q;
}

// 定义一个函数来求解燃料电池一维模型，并打印结果
void solve() {

    // 创建动态数组来存储电势、浓度、电流密度
    double* phi = new double[N]; // 电势
    double* c = new double[N]; // 浓度
    double* i = new double[N]; // 电流密度

    // 给定初始值
    for (int k = 0; k < N; k++) {
        phi[k] = 0.0;
        c[k] = 0.0;
        i[k] = 0.0;
    }

    // 设置边界条件
    phi[1] = phi_a;
    phi[N - 1] = phi_c;

    // 创建动态数组来存储线性方程组的系数矩阵和右端向量
    double* a = new double[N]; // 下对角线元素
    double* b = new double[N]; // 主对角线元素
    double* cc = new double[N]; // 上对角线元素
    double* d = new double[N]; // 右端向量

    // 给定初始值
    for (int k = 0; k < N; k++) {
        a[k] = 0.0;
        b[k] = 0.0;
        cc[k] = 0.0;
        d[k] = 0.0;
    }

    // 定义一个变量来存储误差
    double error = 1e-6;

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

    // 定义一个变量来存储最大迭代次数
    int max_iter = 100;

    // 使用牛顿法迭代求解非线性方程组，直到收敛或达到最大迭代次数
    while (error > 1e-6 && iter < max_iter) {
        // 计算每个节点处的浓度和电流密度
        for (int k = 1; k < N - 1; k++) {
            c[k] = concentration(k, phi[k], c[k]);
            i[k] = current(k, phi[k], c[k]);
        }

        // 根据有限差分法建立线性方程组 Ax = b，其中 x 是电势，A 是系数矩阵，b 是右端向量
        for (int k = 2; k < N - 1; k++) {
            a[k] = -conductivity(k) / dx(k);
            b[k] = conductivity(k) / dx(k) + conductivity(k - 1) / dx(k - 1);
            cc[k] = -conductivity(k - 1) / dx(k - 1);
            d[k] = i[k - 1] - i[k];
        }

        // 处理边界条件，即电势已知的节点
        b[1] = 1.0;
        d[1] = phi_a;
        b[N - 1] = 1.0;
        d[N - 1] = phi_c;

        // 使用追赶法解线性方程组，得到新的电势
        thomas(a, b, cc, d, phi, N);

        // 计算误差，即电流密度在两端的差值的绝对值
        error = abs(i[1] - i[N - 1]);

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

    // 打印结果
    cout << "The simulation results are:" << endl;
    cout << "The iteration number is: " << iter << endl;
    cout << "The error is: " << error << endl;
    cout << "The current density is: " << i[1] << endl;
    cout << "The potential distribution is:" << endl;
    for (int k = 1; k < N; k++) {
        cout << phi[k] << " ";
    }
    cout << endl;

    // 释放动态分配的内存空间
    delete[] phi;
    delete[] c;

    delete[] i;
    delete[] a;
    delete[] b;
    delete[] cc;
    delete[] d;
}

// 定义主函数
int main() {
    // 调用 solve 函数
    solve();
    // 返回 0
    return 0;
}