CFD算例代码：方柱绕流-卡门涡街（C++）

猿核试Bug愁

已于 2025-01-06 15:21:48 修改

阅读量1.1k

点赞数 22

分类专栏： CFD算例代码文章标签： 1024程序员节算法 c++ 数学建模

于 2024-10-24 16:38:35 首次发布

本文链接：https://blog.csdn.net/L_peanut/article/details/143207649

版权

CFD算例代码专栏收录该内容

2 篇文章

订阅专栏

一、物理模型

卡门涡街是经典的计算流体力学（CFD）问题，也是验证流体力学问题计算方法的基准案例之一。方柱绕流是模拟卡门涡街的经典模型之一，主要模拟流体流动过程中绕过方向障碍后的流动变化情况，模型如图1所示，来流方向为水平方向。

二、数学模型

控制方程为：

①连续性方程： $\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$

②x方向动量方程： $\rho(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y})=-\frac{\partial P}{\partial x}+\mu (\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})$

③y方向动量方程： $\rho(\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y})=-\frac{\partial P}{\partial y}+\mu (\frac{\partial^{2}v}{\partial x^{2}}+\frac{\partial^{2}v}{\partial y^{2}})$

其中，u为x方向的速度，v为y方向的速度，P为压力，ρ为流体密度，μ为动力粘度。

三、代码实现

根据功能类型将工程代码分为个部分：

头文件：CubeCircleFlow.h

主程序：CCF.cpp

边界处理：BoundaryTreat.cpp

高斯-赛德尔迭代：GaussSeidel.cpp

初始化数值：InitialValue.cpp

龙格库塔法：RungeKutta.cpp

拉普拉斯变换：LaplaceTransform.cpp

其它：Other.cpp

下面为各功能具体的C++代码实现。

3.1 头文件

主要声明参数和函数。

#ifndef CUBE_CIRCLE_FLOW_H
#define CUBE_CIRCLE_FLOW_H

#include <stdio.h>
#include <string>
#include <time.h>
#include <iostream>
#include <cmath>

using namespace std;

#define PI 3.14159265359
#define INF 10000000
#define N 520

extern int n, NX, NY, i_k, j_k;
extern double Re, dt, eps, scale, h, miu, U, T_max;
extern double a, b, x_k, y_k;
extern double psi[N][N], omega[N][N], vx[N][N], vy[N][N], v[N][N], r[N][N], e[N][N];
extern double xnew[N][N], k1[N][N], k2[N][N], k3[N][N], k4[N][N], temp[N][N];
extern unsigned char img[54+3*N*N];

bool PositionDetermine(int i, int j);                    //判断位置是否属于方柱
void LaplaceTransform(double x, double eps);             //拉普拉斯变换
void InitialValue();                                     //初始化
void GaussSeidel(double x, double p[N][N], int itera, double eps);   //高斯-赛德尔迭代
void RungeKutta(double xnew[N][N], double p[N][N]);      //龙格库塔实现
void BoundaryTreat();                                    //边界处理
void ColorChange(string s, double p[N][N]);              //绘制压力图
void HSV2RGB(double color, double x, double val, double &red, double &green, double &blue);                                                  //RBG转换
void WritePlt(string s);                                 //输出Tecplot文件

#endif

3.2 主程序

程序主体部分。

#include "CubeCircleFlow.h"

int n = 20;                        //边长为1的外接正方形的边的划分数
int NX = 400+1;                    //水平方向网格数
int NY = 200+1;                    //垂直方向网格数
double Re = 100;                   //雷诺数
double dt = 0.001;                 //时间步长
double eps = 1e-4;                 //空间离散计算的精度（收敛标准）
double scale = 10.0;               //颜色转换缩放因子
double h = 1.0/n;                  //空间步长
double miu = 0.001;                //运动粘度
double U = Re*miu;                 //速度
double T_max = 80.0/U;             //问题的特征时间
double a = (double)((double)NX/n); //矩形的长度（将这条边分成a*n块）
double b = (double)((double)NY/n); //矩形的宽度（将这条边分成b*n块）
double x_k = (a-1)/4.0;            //移动块左下角的横坐标
double y_k = (b-1)/2.0;            //移动块左下角的纵坐标
int i_k = (int)(x_k*(double)n+1);  //正方形左下角对应的x方向索引
int j_k = (int)(j_k*(double)n+1);  //正方形左下角对应的y方向索引
double psi[N][N], omega[N][N], vx[N][N], vy[N][N], v[N][N], r[N][N], e[N][N];
double xnew[N][N], k1[N][N], temp[N][N];
unsigned char img[54+3*N*N];

int main(int argc, char *argv[])
{
    int i, j;
    double start, end;

    for(i = 0; i <= NX; i++)
    {
        for(j = 0; j <= NY; j++)
        {
            psi[i][j] = 0;
            omega[i][j] = -1;
        }
    }

    start = clock();
    LaplaceTransform(-0.5, eps);
    end = clock();
    printf("The iteration time is %.3lfs\n", (end-start)/CLOCKS_PER_SEC);
    
    //初次laplace变换后
    ColorChange("./laplace.bmp", psi);

    //初始化流场
    InitialValue();
    ColorChange("./initial.bmp", psi);

    //高斯-赛德尔迭代
    GaussSeidel(-0.5, omega, 200, eps);

    //写出Plt文件
    WritePlt("result.plt");

    return 0;
}

3.3 边界处理

用于流场边界处理。

#include "CubeCircleFlow.h"

void BoundaryTreat()
{
    int i, j;

    for(i = 1; i <= NX; i++)
    {
        vx[i][1] =  0;
        vy[i][1] =  0;
        vx[i][NY] = 0;
        vy[i][NY] = 0;
    }

    for(i = 2; i <= NY-1; i++)
    {
        vx[1][i] =  U;
        vy[1][i] =  0;
        vx[NX][i] = U;
        vy[NX][i] = 0;
    }

    for(i = 2; i <= NX-1; i++)
    {
        for(j = 2; j <= NY-1; j++)
        {
            if(PositionDetermine(i, j))
            {
                vx[i][j] = (psi[i][j+1]-psi[i][j-1])/(2*h);
                vy[i][j] = -(psi[i+1][j]-psi[i-1][j])/(2*h);
            }
            else
            {
                vx[i][j] = 0;
                vy[i][j] = 0;
            }
        }
    }

    for(i = 1; i <= NX; i++)
    {
        for(j = 1;j <= NY; j++)
        {
            v[i][j] = sqrt(vx[i][j]*vx[i][j]+vy[i][j]*vy[i][j]);
        }
    }
}

3.4 高斯-赛德尔迭代

用于迭代更新速度、压力场。

#include "CubeCircleFlow.h"

void GaussSeidel(double x, double p[N][N], int itera, double eps)
{
    int i, j, k, k_max;
    double T, start, end;
    char str[100];

    start = clock();
    k_max = (int)(round(T_max/dt));

    for(k = 0, T = 0; k <= k_max; k++, T += dt)
    {
        for(i = 1; i <= NX; i++)
        {
            omega[i][1]  = 2/(h*h)*(psi[i][1]-psi[i][2]);
            omega[i][NY] = 2/(h*h)*(psi[i][NY]-psi[i][NY-1]);
        }

        for(i = 2; i <= NY-1; i++)
        {  
            omega[1][i]  = 2/(h*h)*(psi[1][i]-psi[2][i]);
            omega[NX][i] = 2/(h*h)*(-psi[NX][i]+psi[NX-1][i]);
        }

        for(i = 0; i < n; i++)
        {
            omega[i_k][j_k+i]     = 2/(h*h)*(psi[i_k][j_k+i]-psi[i_k-1][j_k+i]);
            omega[i_k+n-1][j_k+i] = 2/(h*h)*(psi[i_k+n-1][j_k+i]-psi[i_k+n][j_k+i]);
            omega[i_k+i][j_k]     = 2/(h*h)*(psi[i_k+i][j_k]-psi[i_k+i][j_k-1]);
            omega[i_k+i][j_k+n-1] = 2/(h*h)*(psi[i_k+i][j_k+n-1]-psi[i_k+i][j_k+n]);
        }

        RungeKutta(k1, omega);

        for(i = 2; i <= NX-1; i++)
        {
            for(j = 2; j <= NY-1; j++)
            {
                if(PositionDetermine(i, j))
                {
                    omega[i][j] = omega[i][j]+k1[i][j];
                }
            }
        }

        LaplaceTransform(x, eps);

        if(k%(k_max/itera) == 0)
        {
            printf("k = %d, U = %.2lf, T = %lf\n",k, U, T);
            sprintf(str,"./temp/t = %.3lf.bmp", T);

            end = clock();

            printf("The iteration time is %.3lf\n",(end-start)/CLOCKS_PER_SEC);
            start = clock();

            BoundaryTreat();

            ColorChange(str, p);
        }
    }
}

3.5 初始化数值

用于计算初始化。

#include "CubeCircleFlow.h"

void InitialValue()
{
    int i, j;

    for(i = 1; i <= NX; i++)
    {
        for(j = 1; j <= NY; j++)
        {
                psi[i][j]   = 0;
                omega[i][j] = 0;
        }
    }

    for(i = 1; i <= NX; i++)
    {
        psi[i][NY] = U*b;
    }
    for(i = 1; i <= NY; i++)
    {
        psi[1][i] = psi[NX][i]=U*b*(double)(i-1)/(NY-1);
    }
    for(i = 0; i < n; i++)
    {
        psi[i_k][j_k+i]     = U*b/2;
        psi[i_k+n-1][j_k+i] = U*b/2;
        psi[i_k+i][j_k]     = U*b/2;
        psi[i_k+i][j_k+n-1] = U*b/2;
    }
}

3.6 龙格库塔法

龙格库塔法预测流场速度。

#include "CubeCircleFlow.h"

void RungeKutta(double xnew[N][N], double x[N][N])
{
    int i, j;
    double u_p;

    for(i = 2; i <= NX-1; i++)
    {
        for(j = 2; j <= NY-1; j++)
        {
            if(PositionDetermine(i, j))
            {
                u_p = -(psi[i][j+1]-psi[i][j-1])*(x[i+1][j]-x[i-1][j])/4.0 + (psi[i+1][j]-psi[i-1][j])*(x[i][j+1]-x[i][j-1])/4.0 + miu*(x[i+1][j]+x[i-1][j]+x[i][j+1]+x[i][j-1]-4*x[i][j]);
                xnew[i][j] = dt*u_p/(h*h);
            }
        }
    }
}

3.7 拉普拉斯变换

用于拉普拉斯变换。

#include "CubeCircleFlow.h"

void LaplaceTransform(double x, double eps)
{
    int i, j, k;
    double u_p, sum_k, psi_old;

    sum_k = INF;
    k = 0;

    while(sum_k > eps)
    {
        sum_k = 0.0;
        for(j = 2; j <= NY-1; j++)
        {
            for(i = 2; i <= NX-1; i++)
            {
                if(PositionDetermine(i, j))
                {
                    u_p = (psi[i-1][j]+psi[i+1][j]+psi[i][j-1]+psi[i][j+1]+h*h*omega[i][j])/4.0;
                    psi_old = psi[i][j];
                    psi[i][j] = x*psi[i][j]+(1-x)*u_p;
                    sum_k = sum_k+(psi[i][j]-psi_old)*(psi[i][j]-psi_old);
                }
            }
        }

        sum_k = sqrt(sum_k);
        k++;
    }
}

3.8 其它

用于位置判断、颜色转换、压力云图输出、plt文件写出。

#include "CubeCircleFlow.h"

//判断是否为方形柱
bool PositionDetermine(int i, int j)
{
    if(i_k <= i && i <= i_k+n-1 && j_k <= j && j <= j_k+n-1)
        return 0;
    else
        return 1;
}

void ColorChange(string s, double p[N][N])
{
    int i, j, w, h, x, y;
    int filesize;
    double r, g, b, color;
    double mn = INF, mx = -INF;
    
    w = NX;
    h = NY;

    FILE *f;
    filesize = 54+3*w*h;

    for(i = 1; i <= NX; i++)
    {
        for(j = 1; j <= NY; j++)
        {
            mn = min(mn, p[i][j]);
            mx = max(mx, p[i][j]);
        }
    }

    printf("min = %lf\nmax = %lf\n\n", mn, mx);

    mn = mn/scale;
    mx = mx/scale;

    for(j = 1; j <= NY; j++)
    {
        for(i = 1; i <= NX; i++)
        {
            color = p[i][j];
            color = max(color, mn);
            color = min(color, mx);

            color = 360-360*(color-mn)/(mx-mn);
            color = color*0.8;

            while(color > 360)
                color = color-360;
            while(color < 0)
                color = color+360;

            if(PositionDetermine(i, j))
                HSV2RGB(color, 1, 1, r, g, b);
            else
                HSV2RGB(color, 1, 0, r, g, b);

            r = r*255;
            g = g*255;
            b = b*255;

            x = i-1;
            y = (h-1)-j+1;
            img[(x+y*w)*3+2] = (unsigned char)(r);
            img[(x+y*w)*3+1] = (unsigned char)(g);
            img[(x+y*w)*3+0] = (unsigned char)(b);
        }
    }

    unsigned char bmpfileheader[14] = { 'B','M', 0,0,0,0, 0,0, 0,0, 54,0,0,0 };
    unsigned char bmpinfoheader[40] = {40,0,0,0, 0,0,0,0, 0,0,0,0, 1,0, 24,0 };
    unsigned char bmppad[3] = { 0,0,0 };

    bmpfileheader[2] = (unsigned char)(filesize);
    bmpfileheader[3] = (unsigned char)(filesize >> 8);
    bmpfileheader[4] = (unsigned char)(filesize >> 16);
    bmpfileheader[5] = (unsigned char)(filesize >> 24);

    bmpinfoheader[4] = (unsigned char)(w);
    bmpinfoheader[5] = (unsigned char)(w >> 8);
    bmpinfoheader[6] = (unsigned char)(w >> 16);
    bmpinfoheader[7] = (unsigned char)(w >> 24);
    bmpinfoheader[8] = (unsigned char)(h);
    bmpinfoheader[9] = (unsigned char)(h >> 8);
    bmpinfoheader[10] = (unsigned char)(h >> 16);
    bmpinfoheader[11] = (unsigned char)(h >> 24);

    f = fopen(s.c_str(), "wb");
    fwrite(bmpfileheader, 1, 14, f);
    fwrite(bmpinfoheader, 1, 40, f);
    for(i = 0; i < h; i++)
    {
        fwrite(img+(w*(h-i-1)*3), 3, w, f);
        fwrite(bmppad, 1, (4-(w*3)%4)%4, f);
    }
    fclose(f);
}

void HSV2RGB(double color, double x, double val, double &red, double &green, double &blue)
{
    double i, f, p, q, t;

    if(val == 0)
    {
        red   = 0;
        green = 0;
        blue  = 0;
    }
    else
    {
        color = color/60;
        i = floor(color);
        f = color-i;
        p = val*(1-x);
        q = val*(1-(x*f));
        t = val*(1-(x*(1-f)));
        if(i == 0)      {red=val; green=t; blue=p;}
        else if(i == 1) {red=q; green=val; blue=p;}
        else if(i == 2) {red=p; green=val; blue=t;}
        else if(i == 3) {red=p; green=q; blue=val;}
        else if(i == 4) {red=t; green=p; blue=val;}
        else if(i == 5) {red=val; green=p; blue=q;}
    }
}

void WritePlt(string s)
{
    int i, j;
    FILE *f=fopen(s.c_str(), "w");

    for(i = 1; i <= NX; i++)
    {
        for(j = 1; j <= NY; j++)
        {
            fprintf(f, "%lf", psi[i][j]);
        }
        putc(10, f);
    }

    for(i = 1; i <= NX; i++)
    {
        for(j = 1; j <= NY; j++)
        {
            fprintf(f, "%lf", omega[i][j]);
        }
        putc(10, f);
    }
    
    fclose(f);
}

四、测试结果

为方便代码编译，编写Makefile文件，使用make命令进行编译。

SRC=CCF.cpp BoundaryTreat.cpp GaussSeidel.cpp InitialValue.cpp LaplaceTransform.cpp RungeKutta.cpp Other.cpp
OBJ=$(SRC:.cpp=.o)
TARGET=CCF

CC=g++
CFLAGS=-O2

$(TARGET):$(OBJ)
	$(CC) $(CFLAGS) -o $@ $^
%.o:%.cpp
	$(CC) $(CFLAGS) -c $? -o $@

.PHONY:clean

clean:
	rm -rf $(TARGET) $(OBJ)

编译：