二维结构网格单相流SIMPLE 算法的C语言程序代码

以下是一个简单的二维结构网格单相流SIMPLE算法的C语言程序代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NX 21 // X方向网格数
#define NY 21 // Y方向网格数
#define L 1.0 // 计算区域长度
#define H 1.0 // 计算区域高度
#define DX (L / (NX - 1)) // X方向网格大小
#define DY (H / (NY - 1)) // Y方向网格大小
#define U 1.0 // 初始速度
#define V 0.0 // 初始速度
#define RHO 1.0 // 流体密度
#define MU 0.1 // 流体动力粘度
#define RE (U * H / MU) //雷诺数
#define DT (0.01 * DX * DX / MU) // 时间步长
#define MAX_STEP 10000 // 最大迭代次数
#define TOL 1e-6 // 迭代收敛精度

double u[NX][NY], v[NX][NY]; // 速度场
double p[NX][NY]; // 压力场
double b[NX][NY]; // Poisson方程右端项

// 边界条件处理
void boundary_conditions() {
    // 上下边界
    for (int i = 0; i < NX; i++) {
        u[i][0] = U;
        v[i][0] = 0.0;
        p[i][0] = p[i][1];
        u[i][NY - 1] = 0.0;
        v[i][NY - 1] = 0.0;
        p[i][NY - 1] = p[i][NY - 2];
    }
    // 左右边界
    for (int j = 0; j < NY; j++) {
        u[0][j] = 0.0;
        v[0][j] = 0.0;
        p[0][j] = p[1][j];
        u[NX - 1][j] = 0.0;
        v[NX - 1][j] = 0.0;
        p[NX - 1][j] = p[NX - 2][j];
    }
}

// Poisson方程求解
void pressure_poisson() {
    double delta;
    int step = 0;
    while (step < MAX_STEP) {
        delta = 0.0;
        for (int i = 1; i < NX - 1; i++) {
            for (int j = 1; j < NY - 1; j++) {
                double t = ((p[i + 1][j] + p[i - 1][j]) * DY * DY +
                            (p[i][j + 1] + p[i][j - 1]) * DX * DX -
                            b[i][j] * DX * DX * DY * DY) /
                           (2 * (DX * DX + DY * DY));
                delta += fabs(t - p[i][j]);
                p[i][j] = t;
            }
        }
        if (delta < TOL)
            break;
        step++;
    }
}

// 速度场更新
void velocity_update() {
    for (int i = 1; i < NX - 1; i++) {
        for (int j = 1; j < NY - 1; j++) {
            u[i][j] = u[i][j] - DT / RHO * (p[i + 1][j] - p[i - 1][j]) / (2 * DX);
            v[i][j] = v[i][j] - DT / RHO * (p[i][j + 1] - p[i][j - 1]) / (2 * DY);
        }
    }
}

// Poisson方程右端项计算
void calculate_b() {
    for (int i = 1; i < NX - 1; i++) {
        for (int j = 1; j < NY - 1; j++) {
            b[i][j] = RHO * (1.0 / DT *
                             ((u[i + 1][j] - u[i - 1][j]) / (2 * DX) +
                              (v[i][j + 1] - v[i][j - 1]) / (2 * DY)) -
                             ((u[i + 1][j] - u[i - 1][j]) / (2 * DX)) *
                             ((u[i + 1][j] - u[i - 1][j]) / (2 * DX)) -
                             2 * ((u[i][j + 1] - u[i][j - 1]) / (2 * DY) *
                                  (v[i + 1][j] - v[i - 1][j]) / (2 * DX)) -
                             ((v[i][j + 1] - v[i][j - 1]) / (2 * DY)) *
                             ((v[i][j + 1] - v[i][j - 1]) / (2 * DY)));
        }
    }
}

// 输出结果
void output_results(int step) {
    char filename[100];
    sprintf(filename, "result_%04d.dat", step);
    FILE *fp = fopen(filename, "w");
    for (int i = 0; i < NX; i++) {
        for (int j = 0; j < NY; j++) {
            fprintf(fp, "%.2f %.2f %.5f %.5f %.5f\n", i * DX, j * DY, u[i][j], v[i][j], p[i][j]);
        }
        fprintf(fp, "\n");
    }
    fclose(fp);
}

int main() {
    // 初始化
    for (int i = 0; i < NX; i++) {
        for (int j = 0; j < NY; j++) {
            u[i][j] = U;
            v[i][j] = V;
            p[i][j] = 0.0;
        }
    }
    // 边界条件处理
    boundary_conditions();
    // 迭代求解
    for (int step = 0; step < MAX_STEP; step++) {
        // Poisson方程右端项计算
        calculate_b();
        // Poisson方程求解
        pressure_poisson();
        // 速度场更新
        velocity_update();
        // 边界条件处理
        boundary_conditions();
        // 输出结果
        output_results(step);
    }
    return 0;
}

以上代码中的SIMPLE算法使用了标准的高斯-赛德尔迭代方法求解Poisson方程，具体实现可以参考代码中的pressure_poisson()函数。输出结果的函数output_results()将结果保存在文件中，可以用Paraview等软件进行后处理和可视化。