CFD算例代码：顶盖驱动方腔流（C++）

猿核试Bug愁

已于 2024-10-23 15:29:26 修改

阅读量3.7k

点赞数 28

CC 4.0 BY-SA版权

分类专栏： # CFD算例代码文章标签：算法数学建模 c++

于 2024-10-23 15:24:21 首次发布

本文链接：https://blog.youkuaiyun.com/L_peanut/article/details/143160147

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

2 篇文章

订阅专栏

4.1 顶盖速度U=0.05，Re=500

4.2 顶盖速度U=0.3，Re=3000

4.3 顶盖速度U=1.0，Re=10000

一、物理模型

顶盖驱动方腔流（Lid-driven Cavity Flow）是经典的计算流体力学（CFD）问题，也是验证流体力学问题计算方法的基准案例之一。二维顶盖驱动方腔流问题的物理模型由一个二维方形空腔构成，其上壁面具有水平方向的速度，其它三个壁面为静止壁面（速度为0），模型如图1所示。流体力学实验表明，该形式的运动会产生一种流动现象，特征是：空腔中心会形成一个大涡，空腔角落形成一些小涡，而雷诺数的大小将会影响涡的大小和数量。流动特征如图2所示。

二、数学模型

2.1 流动控制方程

采用不可压缩定常流Navier-Stokes方程：

$\rho (\mathbf{u}\cdot \triangledown )\mathbf{u}=- \triangledown p+\mu \triangledown^{2}\mathbf{u}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(1)$

将公式（1）中速度、压力、长度参数无量纲化后得到不可压缩定常流Navier-Stokes方程：

$(\frac{u}{v}\cdot L\triangledown )\frac{u}{v}=-\triangledown \frac{p}{\rho v^{2}}+\frac{1}{Re}(L\triangledown )^{2}\frac{u}{v}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(2)$

将公式（2）写成x和y方向的二维方程组，可得：

$\left\{\begin{matrix} \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0\\ \frac{\partial u}{\partial t}+\frac{\partial (uu)}{\partial x}+\frac{\partial (uv)}{\partial y}=-\frac{\partial p}{\partial x}+\frac{1}{Re}(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})\\ \frac{\partial v}{\partial t}+\frac{\partial (vu)}{\partial x}+\frac{\partial (vv)}{\partial y}=-\frac{\partial p}{\partial y}+\frac{1}{Re}(\frac{\partial^{2}v}{\partial x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}) \end{matrix}\right.\;\;\;(3)$

2.2 离散方法

1、求解域离散

计算域为方形，边长 $L=0.1m\;(0\leqslant x\leqslant 0.1,\;0\leqslant y\leqslant 0.1)$ ，采用二维均匀网格划分，将计算域离散为 $n\times n$ 的均匀网格，空间步长为 $\Delta x=\Delta y =\Delta h = L/n=0.1/n$ ，离散后的网格模型如下图所示（n=200）。

2、控制方程离散

模型的边界条件为：

$U|_{x=0}=U|_{x=L}=U|_{y=0}=0,\;\;\;\;\;\;\;U|_{y=L}=1$

$V|_{x=0}=V|_{x=L}=V|_{y=0}=V|_{y=L}=0$

$\frac{\partial P}{\partial x}|_{x=0}=\frac{\partial P}{\partial x}|_{x=L}=\frac{\partial P}{\partial y}|_{y=0}=\frac{\partial P}{\partial y}|_{y=L}=0$

使用交错网格，将U和V的边界条件视为虚拟节点与内部一层节点之和的1/2等于0，即虚拟节点与内部节点互为相反数。P的边界条件在网格中表现为第一层和第二层（或最后一层和倒数第二层）相等。

差分方程：

$u^{n^{*}}_{i+1/2,j}=u^{n}_{i+1/2,j}+A\Delta t-\frac{\Delta t}{\Delta x}(p^{n}_{i+1,j}-p^{n}_{i,j})\;\;\;\;\;\;\;(4a)$

$v^{n^{*}}_{i,j+1/2}=v^{n}_{i,j+1/2}+B\Delta t-\frac{\Delta t}{\Delta y}(p^{n}_{i,j+1}-p^{n}_{i,j})\;\;\;\;\;\;\;(4b)$

将假设的U和V代入差分方程（4a）和（4b），得到当前的U*和V*，如下式：

$A=-[\frac{(u^{2})^{n}_{i+3/2,j}-(u^{2})^{n}_{i-1/2,j}}{2\Delta x}+\frac{u^{n}_{i+1/2,j+1}\bar{v}-u^{n}_{i+1/2,j-1}v}{2\Delta y}]+\frac{1}{Re}[\frac{u^{n}_{i+3/2,j}-2u^{n}_{i+1/2,j}+u^{n}_{i-1/2,j}}{(\Delta x)^{2}}+\frac{u^{n}_{i+1/2,j+1}-2u^{n}_{i+1/2,j}+u^{n}_{i+1/2,j-1}}{(\Delta y)^{2}}]\;\;\;\;\;\;\;(5a)$

$B=-[\frac{v^{n}_{i+1,j+1/2}\bar{u}-v^{n}_{i-1,j+1/2}u}{2\Delta x}+\frac{(v^{2})^{n}_{i+1/2,j+3/2}-(v^{2})^{n}_{i,j+3/2}}{2\Delta y}]+\frac{1}{Re}[\frac{v^{n}_{i+1,j+1/2}-2v^{n}_{i,j+1/2}+v^{n}_{i-1,j+1/2}}{(\Delta x)^{2}}+\frac{v^{n}_{i,j+3/2}-2v^{n}_{i,j+1/2}+v^{n}_{i,j-1/2}}{(\Delta y)^{2}}]\;\;\;\;\;\;\;(5b)$

其中：

$\bar{v}=\frac{1}{2}(v_{i,j+1/2}+v_{i+1,j+1/2}),\;\;v=\frac{1}{2}(v_{i,j-1/2}+v_{i+1,j-1/2})$

$u=\frac{1}{2}(u_{i-1/2,j}+u_{i-1/2,j+1}),\;\;\bar{u}=\frac{1}{2}(u_{i+1/2,j}+u_{i+1/2,j+1})$

2.3 迭代方法

迭代方法分别采用雅各比（Jacobi）迭代和高斯-赛德尔（Gauss-Seidel）迭代来求解，这两种迭代方法的具体理论推导请参考博主的另一篇博文：↓传送门↓

偏微分方程算法之迭代法（番外篇）_迭代法解偏微分方程-优快云博客https://blog.youkuaiyun.com/L_peanut/article/details/138302224 将U*和V*代入公式（6），可得到P'，利用雅各比迭代或高斯-赛德尔迭代求解方程组：

$P^{'}_{i,j}=-\frac{1}{a}(bP^{'}_{i+1,j}+bP^{'}_{i-1,j}+cP^{'}_{i,j+1}+cP^{'}_{i,j-1}+d)\;\;\;\;\;\;\;(6)$

其中：

$a=2[\frac{\Delta t}{(\Delta x)^{2}}+\frac{\Delta t}{(\Delta y)^{2}}]$

$b=-\frac{\Delta t}{(\Delta x)^{2}}$

$c=-\frac{\Delta t}{(\Delta y)^{2}}$

$d=\frac{1}{\Delta x}[(u^{*})_{i+1/2,j}-(u^{*})_{i-1/2,j}]+\frac{1}{\Delta y}[(v^{*})_{i,j+1/2}-(v^{*})_{i,j-1/2}]$

将P'代入 $P=P^{*}+0.1\times P^{'}$ 中，得到当前时刻准确的P值，再将P'代入公式（7a）和（7b），得到当前时刻准确的U和V值。

$U^{n}_{i+1/2,j}=U^{n^{*}}_{i+1/2,j}-\frac{\Delta t}{\Delta x}(P^{'}_{i+1,j}-P^{'}_{i,j})\;\;\;\;\;\;\;(7a)$

$V^{n}_{i,j+1/2}=V^{n^{*}}_{i,j+1/2}-\frac{\Delta t}{\Delta y}(P^{'}_{i,j+1}-P^{'}_{i,j})\;\;\;\;\;\;\;(7b)$

再将当前时刻准确的U、V和修正的P代入差分方程（4a）和（4b）就可以得到下一个时刻的U和V，如此循环。

三、代码实现

为方便测试采用不同迭代方法的算法，按照功能将工程代码划分为五部分：

头文件：CavityFlow.h

主程序：LidDrivenCavityFlow.cpp

雅各比迭代：Jacobi.cpp

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

结果输出：OutputPlt.cpp

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

3.1 头文件

主要声明参数和函数。

#ifndef CAVITY_FLOW
#define CAVITY_FLOW

#include <iostream>
#include <fstream>
#include <string>
#include <cmath>

using namespace std;

#define L 0.1                        //方腔边长，单位：米(m)
#define N 2000                       //外场大小
#define BD 201                       //计算域边界
#define T 0.00005                    //时间步长
#define UTOP 0.3                     //顶盖运动速度，单位：米/秒(m/s)
#define RXFACTOR 1.0                 //松弛因子，等于1.0时为高斯-赛德尔迭代，大于1.0 
                                     //时为超松弛迭代，小于1.0时为亚松弛迭代
#define VISCOSITY 0.00001            //粘性系数，单位:帕斯卡·秒(Pa·s)
#define EPSILON1 0.00001             //全局收敛标准
#define EPSILON2 0.000001            //松弛迭代收敛标准
#define STEP 10000                    //计算域整体网格迭代上限
#define SINGLESTEP 1000               //计算域单个网格迭代上限

ofstream WriteResidual(double Re);
void WriteResult(double **U, double **V, double **P, double Re);
#ifdef GS
double GaussSeidel(double **U1, double **V1, double **P2, double **P3, int &itera1);
#else
double Jacobi(double **U1, double **V1, double **P2, double **P3, int &itera1);
#endif

#endif

3.2 主程序

程序主体部分。

/*---------------------------------------------------------------
				Lid-Driven Cavity Flow

						U ——→
	    ||---------------------------------|| Boundary-condition:
	    ||	↑							   ||     	U(x=0)=0
	    ||	| y	   					       ||	  	U(x=L)=0
	    ||								   ||	  	U(y=0)=0
	    ||								   ||	  	U(y=L)=1
	    ||								   ||	  	V(x=0)=0
	    ||								   ||	  	V(x=L)=0
	    ||								   ||	  	V(y=0)=0
	    ||								   ||	  	V(y=L)=0
	    ||	     x   					   ||	  
	    ||	     ——→ 					   ||	  	αP/αx =0
		||=================================||	  	αP/αy =0
		   ←————————————— L —————————————→ 

	  	Two Dimension Incompressible N-S equations

  αu/αt+u*αu/αx+v*αu/αy = -1/ρ*αp/αx+μ*α2u/(αx2)+μ*α2u/(αy2)
  αv/αt+u*αv/αx+v*αv/αy = -1/ρ*αp/αy+μ*α2v/(αx2)+μ*α2v/(αy2)
	     αu/αx + αv/αy  = 0

  "u","v" 分别为x,y方向上的速度;
  "t"为时间;
  "ρ"为密度,不可压缩流体的密度为常数;
  "μ"为流体粘性系数;
  "p"为压强.
---------------------------------------------------------------*/
#include "CavityFlow.h"

int main(int argc, char *argv[])
{
    int i, j, n, m;
    int itera1, itera2;
    double x, y;
    double ux, uy, vx, vy;
    double a, b, c, d, e;
    double Re, k2, sum2;
    double **U, **V, **U1, **V1, **U2, **V2;
    double **P, **P1, **P2, **P3;

    //分配空间
    U = (double **)malloc(sizeof(double*)*(N+1));
    U1 = (double **)malloc(sizeof(double*)*(N+1));
    U2 = (double **)malloc(sizeof(double*)*(N+1));
    V = (double **)malloc(sizeof(double*)*(N+1));
    V1 = (double **)malloc(sizeof(double*)*(N+1));
    V2 = (double **)malloc(sizeof(double*)*(N+1));
    P = (double **)malloc(sizeof(double*)*(N+1));
    P1 = (double **)malloc(sizeof(double*)*(N+1));
    P2 = (double **)malloc(sizeof(double*)*(N+1));
    P3 = (double **)malloc(sizeof(double*)*(N+1));

    for(i = 0; i <= N; i++)
    {
        U[i] = (double *)malloc(sizeof(double)*(N+1));
        U1[i] = (double *)malloc(sizeof(double)*(N+1));
        U2[i] = (double *)malloc(sizeof(double)*(N+1));
        V[i] = (double *)malloc(sizeof(double)*(N+1));
        V1[i] = (double *)malloc(sizeof(double)*(N+1));
        V2[i] = (double *)malloc(sizeof(double)*(N+1));
        P[i] = (double *)malloc(sizeof(double)*(N+1));
        P1[i] = (double *)malloc(sizeof(double)*(N+1));
        P2[i] = (double *)malloc(sizeof(double)*(N+1));
        P3[i] = (double *)malloc(sizeof(double)*(N+1));
    }

    //x、y为空间步长，Re为雷诺数
    x = L/(BD-1), y = L/(BD-1);
    a = (T/(x*x) + T/(y*y))*2;
    b = -T/(x*x);
    c = -T/(y*y);
    Re = (1.0*UTOP*L)/VISCOSITY;
    itera1 = 0, itera2 = 0, sum2 = 0.0;
    cout<<"a="<<a<<", b="<<b<<", c="<<c<<", Re="<<Re<<endl;

#ifdef GS
    cout<<"Iteration Function:   GaussSeidel"<<endl;
#else
    cout<<"Iteration Function:   Jacobi"<<endl;
#endif

    //对假设的压力场P*进行初始化
    for(i = 1; i <= BD; i++)
    {
        for(j = 1; j <= BD; j++)
        {
            P[i][j] = 0;
        }
    }
    //求解域内u*,v*的初值
    for(i = 1; i <=BD-1; i++)
    {
        for(j = 1; j <= BD-1; j++)
        {
            U[i][j] = 0;
            V[i][j] = 0;
        }
    }
    //边界上u*,v*的初值
    for(i = 0; i <= BD; i++)
    {
        U[i][1] =  0;
        U[i][BD] = UTOP;
        V[1][i] =  0;
        V[BD][i] = 0;
    }
    for(i = 1; i <= BD; i++)
    {
        U[0][i] =  0 - U[1][i];
        U[BD][i] = 0 - U[BD-1][i];
        V[i][0] =  0 - V[i][1];
        V[i][BD] = 0 - V[i][BD-1];
    }
    for(i = 0; i < BD; i++)
    {
        U[i][0] =    0 - U[i][2];
        U[i][BD+1] = 2 - U[i][BD-1];
        V[0][i] =    0 - V[2][i];
        V[BD+1][i] = 0 - V[BD-1][i];
    }

    //残差记录
    ofstream fout = WriteResidual(Re);

    /*----------------------------------------------------------------------------
	计算预测流场和压力场:
	1、U(n+1)[i][j] = U(n)[i][j] - (U(n)[i+1][j]^2 - U(n)[i][j]^2)/x*t - (ux-uy)/y*t + t/Re*((U(n)[i+1][j]-2*U(n)[i][j]+U(n)[i-1][j])/x^2 + (U(n)[i][j+1]-2*U(n)[i][j]+U(n)[i][j-1])/y^2) - t/x*(P(n)[i+1][j]-P[i][j])
	2、V(n+1)[i][j] = V(n)[i][j] - (V(n)[i+1][j]^2 - V(n)[i][j]^2)/y*t - (vx-vy)/x*t + t/Re*((V(n)[i+1][j]-2*V(n)[i][j]+V(n)[i-1][j])/x^2 + (V(n)[i][j+1]-2*V(n)[i][j]+V(n)[i][j-1])/y^2) - t/y*(P(n)[i+1][j]-P[i][j])
	3、P(n+1)([i][j] = -1/a*(b*P(n)[i+1][j] + b*P(n)[i-1][j] + c*P(n)[i][j+1] + c*P(n)[i][j-1] + 1/x*(U(n+1)[i][j]-U(n+1)[i-1][j]) + 1/y*(V(n+1)[i][j]-V(n+1)[i][j-1]))
	----------------------------------------------------------------------------*/
	for(n = 1; n <= STEP; n++)
    {
        /*-----------------------------------------------------------------------
		-----用假设的压力场求解动量方程，得到下一步的未修正流场U(n+1)*, V(n+1)* ，但不满足连续性方程-----
		------------------------------------------------------------------------*/
        for(i = 1; i <= BD-1; i++)
        {
            for(j = 1; j <= BD; j++)
            {
                ux=(U[i][j+1]+U[i][j])*(V[i+1][j]+V[i][j])/4;
				uy=(U[i][j]+U[i][j-1])*(V[i+1][j-1]+V[i][j-1])/4;
   				U1[i][j]=U[i][j]-(U[i+1][j]*U[i+1][j]-U[i][j]*U[i][j])/x*T-(ux-uy)/y*T+1/Re*((U[i+1][j]-2*U[i][j]+U[i-1][j])/(x*x)+(U[i][j+1]-2*U[i][j]+U[i][j-1])/(y*y))*T-T/x*(P[i+1][j]-P[i][j]);
            }
        }
        for(i = 1; i <= BD; i++)
        {
            for(j = 1; j <= BD-1; j++)
            {
                vx=(U[i][j+1]+U[i][j])*(V[i+1][j]+V[i][j])/4;
		    	vy=(U[i-1][j+1]+U[i-1][j])*(V[i-1][j]+V[i][j])/4;
    			V1[i][j]=V[i][j]-(V[i][j+1]*V[i][j+1]-V[i][j]*V[i][j])/y*T-(vx-vy)/x*T+1/Re*((V[i+1][j]-2*V[i][j]+V[i-1][j])/(x*x)+(V[i][j+1]-2*V[i][j]+V[i][j-1])/(y*y))*T-T/y*(P[i][j+1]-P[i][j]);
            }
        }

        //初始化迭代所需的压力场数组
        for(i = 0; i <= BD+1; i++)
        {
            for(j = 0; j <= BD+1; j++)
            {
                P1[i][j] = 0;
                P2[i][j] = 0;
                P3[i][j] = 0;
            }
        }
        d = 1/x*(U1[10][11]-U1[9][11])+1/y*(V1[10][11]-V1[10][10]);
        P1[10][11] = -b/a*P1[11][11]-b/a*P1[9][11]-c/a*P1[10][12]-c/a*P1[10][10]-d/a;

        /*-----------------------------------------------------------------------------------------
		--------------用得到的未修正U(n+1)*、V(n+1)*和P1计算全局压力场P3，但P3不满足动量方程-----------
		------------------此处计算获得的P3将作为迭代求解泊松方程的初值，从而获得修正后的P3--------------
		-----------------------------(此处第一步采用松弛迭代计算P3)----------------------------------
		------------------------------------------------------------------------------------------*/
        for(i = 2; i <= BD-1; i++)
        {
            for(j = 2; j <= BD-1; j++)
            {
                if(i == 2 && j != 2 && j != BD-1)
                {
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
                    P2[i][j] = P1[i][j] - RXFACTOR * (P1[i][j] + b/(a+b)*P1[i+1][j] + c/(a+b)*P1[i][j+1] + c/(a+b)*P2[i][j-1] + d/(a+b));
                }
                else if(i == BD-1 && j != 2 && j != BD-1)
                {
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P1[i][j] - RXFACTOR * (P1[i][j] + b/(a+b)*P2[i-1][j] + c/(a+b)*P1[i][j+1] + c/(a+b)*P2[i][j-1] + d/(a+b));
                }
                else if(j == 2 && i != 2 && i != BD-1)
                {
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P1[i][j] - RXFACTOR * (P1[i][j] + b/(a+c)*P1[i+1][j] + b/(a+c)*P2[i-1][j] + c/(a+c)*P1[i][j+1] + d/(a+c));
                }
                else if(j == BD-1 && i != 2 && i != BD-1)
                {
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P1[i][j] - RXFACTOR * (P1[i][j] + b/(a+c)*P1[i+1][j] + b/(a+c)*P2[i-1][j] + c/(a+c)*P2[i][j-1] + d/(a+c));
                }
                else if(i == 2 && j == 2)
                {
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P1[i][j] - RXFACTOR * (P1[i][j] + b/(a+b+c)*P1[i+1][j] + c/(a+b+c)*P1[i][j+1] + d/(a+b+c));
                }
                else if(i == 2 && j == BD-1)
                {
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P1[i][j] - RXFACTOR * (P1[i][j] + b/(a+b+c)*P1[i+1][j] + c/(a+b+c)*P2[i][j-1] + d/(a+b+c));
                }
                else if(i == BD-1 && j == 2)
                {
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P1[i][j] - RXFACTOR * (P1[i][j] + b/(a+b+c)*P2[i-1][j] - c/(a+b+c)*P1[i][j+1] + d/(a+b+c));
                }
                else if(i == BD-1 && j ==BD-1)
                {
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P1[i][j] - RXFACTOR * (P1[i][j] + b/(a+b+c)*P2[i-1][j] + c/(a+b+c)*P2[i][j-1] + d/(a+b+c));
                }
                else
                {
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P1[i][j] - RXFACTOR * (P1[i][j] + b/a*P1[i+1][j] + b/a*P2[i-1][j] + c/a*P1[i][j+1] + c/a*P2[i][j-1] + d/a);
                }
            }
        }

        for(i = 2; i <= BD-1; i++)
        {
            for(j = 2; j <= BD-1; j++)
            {
                if(i == 2 && j != 2 && j != BD-1)
                {
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b)*P2[i+1][j] + c/(a+b)*P2[i][j+1] + c/(a+b)*P3[i][j-1] + d/(a+b));
                }
                else if(i == BD-1 && j !=2 &&j != BD-1)
				{
					d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b)*P3[i-1][j] + c/(a+b)*P2[i][j+1] + c/(a+b)*P3[i][j-1] + d/(a+b));
				}
				else if(j == 2 && i !=2 && i != BD-1)
				{
					d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+c)*P2[i+1][j] + b/(a+c)*P3[i-1][j] + c/(a+c)*P2[i][j+1] + d/(a+c));
				}
				else if(j == BD-1 && i !=2 && i != BD-1)
				{
					d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+c)*P2[i+1][j] + b/(a+c)*P3[i-1][j] + c/(a+c)*P3[i][j-1] + d/(a+c));
				}
				else if(i == 2 && j == 2)
				{
					d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b+c)*P2[i+1][j] + c/(a+b+c)*P2[i][j+1] + d/(a+b+c));
				}
				else if(i == 2 && j == BD-1)
				{
					d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b+c)*P2[i+1][j] + c/(a+b+c)*P3[i][j-1] + d/(a+b+c));
				}
				else if(i == BD-1 && j == 2)
				{
					d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b+c)*P3[i-1][j] + c/(a+b+c)*P2[i][j+1] + d/(a+b+c));
				}
				else if(i == BD-1 && j == BD-1)
				{
					d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b+c)*P3[i-1][j] + c/(a+b+c)*P3[i][j-1] + d/(a+b+c));
				}
				else
				{
					d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/a*P2[i+1][j] + b/a*P3[i-1][j] + c/a*P2[i][j+1] + c/a*P3[i][j-1] + d/a);
				}
            }
        }
    
        /*-------------------------------------------------------------------------------------------------------------
		------------以未修正流场U、V计算得到的P3作为当前压力，进行高斯-赛德尔/雅可比迭代，直到满足迭代步数或者收敛标准---------
		-------------------------------------------------------------------------------------------------------------*/
		//默认采用雅可比迭代法(可通过-DGS选择使用高斯-赛德尔迭代，收敛速度更快)
        #ifdef GS
        e = GaussSeidel(U1, V1, P2, P3, itera1);
        #else
        e = Jacobi(U1, V1, P2, P3, itera1);
        #endif

        //修正压力场P3边界条件处理
        for(i = 1; i <= BD; i++)
        {
            P3[i][1] = P3[i][2];
            P3[i][BD] = P3[i][BD-1];
        }
        for(i = 2; i <= BD-1; i++)
        {
            P3[1][i] = P3[2][i];
            P3[BD][i] = P3[BD-1][i];
        }

        //用修正压力场P3计算速度修正值U',V'
        for(i = 1; i <= BD-1; i++)
        {
            for(j = 1; j <= BD; j++)
            {
                U2[i][j] = -T/x*(P3[i+1][j]-P3[i][j]);
            }
        }
        for(i = 1; i <= BD; i++)
        {
            for(j = 1; j <= BD-1; j++)
            {
                V2[i][j] = -T/y*(P3[i][j+1]-P3[i][j]);
            }
        }

        //准确压力=预测压力+修正压力, P(n+1)=P*+P'
        for(i = 1; i < BD+1; i++)
        {
            for(j = 1; j < BD+1; j++)
            {
                P[i][j] = P[i][j] + P3[i][j];
            }
        }
        //准确速度=预测速度+修正速度, U(n+1)=U*+U'
        for(i = 1; i <= BD-1; i++)
        {
            for(j = 1; j <= BD; j++)
            {
                U[i][j] = U1[i][j] + U2[i][j];
            }
        }
        for(i = 1; i <= BD; i++)
        {
            for(j = 1; j <= BD-1; j++)
            {
                V[i][j] = V1[i][j] + V2[i][j];
            }
        }

        //处理边界速度
        for(i = 0; i <= BD; i++)
        {
            V[1][i] =  0;
            V[BD][i] = 0;
            U[i][1] =  0;
            U[i][BD] = 1;
        }
        for(i = 1; i <= BD; i++)
        {
            U[0][1] =  0 - U[1][i];
            U[BD][i] = 0 - U[BD-1][i];
        }
        for(i = 0; i <= BD; i++)
        {
            U[i][0] =    0 - U[i][2];
            U[i][BD+1] = 2 - U[i][BD-1];
        }
        for(i = 1; i <= BD; i++)
        {
            V[i][0] =  0 - V[i][1];
            V[i][BD] = 0 - V[i][BD-1];
        }
        for(i = 0; i <= BD; i++)
        {
            V[0][i] =    0 - V[2][i];
            V[BD+1][i] = 0 - V[BD-1][i];
        }

        for(i = 2; i <= BD-1; i++)
        {
            for(j = 2; j <= BD-1; j++)
            {
                d = abs(1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]));
				sum2 = sum2+d*d;
            }
        }

        k2 = sqrt(sum2)/((BD-2)*(BD-2));
        if(k2 <= EPSILON1 && n != 1)
        {
            n = STEP + 1;
        }
        sum2 = 0;
        itera2++;

        fout<<n<<" "<<e<<" "<<k2<<"\n";
        cout<<"n="<<n<<", e="<<e<<", k2="<<k2<<endl;
    }
    fout<<endl;
    fout.close();
    cout<<"itera1="<<itera1<<", itera2="<<itera2<<endl;

    //输出网格及结果tecplot文件
    WriteResult(U, V, P, Re);

    //释放内存
    for(i = 0;i <= N;i++)
	{
		free(U[i]);free(U1[i]);free(U2[i]);
		free(V[i]);free(V1[i]);free(V2[i]);
		free(P[i]);free(P1[i]);free(P2[i]);free(P3[i]);
	}
	free(U);free(U1);free(U2);
	free(V);free(V1);free(V2);
	free(P);free(P1);free(P2);free(P3);

    return 0;
}

3.3 雅各比迭代

采用雅各比迭代法更新压力、速度场。

/*---------------------------雅克比迭代-------------------------------
------------P2 = -b/a*P3 - b/a*P3 - c/a*P3 - c/a*P3 - d/a------------
-------------------------------------------------------------------*/

#include "CavityFlow.h"

double Jacobi(double **U1, double **V1, double **P2, double **P3, int &itera1)
{
	int i, j, m;
	double a, b, c, d, e;
    double x, y;
	double k1, sum1;

    x = L/(BD-1), y = L/(BD-1);
	a = (T/(x*x)+T/(y*y))*2;
	b = -T/(x*x);
	c = -T/(y*y);

	for(m = 1; m <= SINGLESTEP; m++)
	{
		for(i = 2; i <= BD-1; i++)
		{
			for(j = 2; j <= BD-1; j++)
			{    
				if(i == 2 && j !=2 && j != BD-1)
			 	{   
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = -b/(a+b)*P3[i+1][j] - c/(a+b)*P3[i][j+1] - c/(a+b)*P3[i][j-1] - d/(a+b);
				}
				else if(i == BD-1 && j !=2 && j != BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = -b/(a+b)*P3[i-1][j] - c/(a+b)*P3[i][j+1] - c/(a+b)*P3[i][j-1] - d/(a+b);
				}
				else if(j == 2 && i !=2 && i != BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = -b/(a+c)*P3[i+1][j] - b/(a+c)*P3[i-1][j] - c/(a+c)*P3[i][j+1] - d/(a+c);
				}
				else if(j == BD-1 && i !=2 && i != BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = -b/(a+c)*P3[i+1][j] - b/(a+c)*P3[i-1][j] - c/(a+c)*P3[i][j-1] - d/(a+c);
				}
				else if(i == 2 && j == 2)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = -b/(a+b+c)*P3[i+1][j] - c/(a+b+c)*P3[i][j+1] - d/(a+b+c);
				}
				else if(i == 2 && j == BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = -b/(a+b+c)*P3[i+1][j] - c/(a+b+c)*P3[i][j-1] - d/(a+b+c);
				}
				else if(i == BD-1 && j == 2)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = -b/(a+b+c)*P3[i-1][j] - c/(a+b+c)*P3[i][j+1] - d/(a+b+c);
				}
				else if(i == BD-1 && j == BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = -b/(a+b+c)*P3[i-1][j] - c/(a+b+c)*P3[i][j-1] - d/(a+b+c);
				}
				else
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = -b/a*P3[i+1][j] - b/a*P3[i-1][j] - c/a*P3[i][j+1] - c/a*P3[i][j-1] - d/a;
				}
			}
		}
		
		for(i = 2; i <= BD-1 ; i++)
		{
			for(j = 2; j <= BD-1; j++)
			{
				k1 = abs(P3[i][j]-P2[i][j]);
				sum1 = sum1+k1*k1;
			}
		}

		e = sqrt(sum1)/((BD-2)*(BD-2));
		if(e <= EPSILON2 && m != 1)
		{
			m = SINGLESTEP+1;
		}
		sum1 = 0;
		itera1++;
	}
	return e;
}

3.4 高斯-赛德尔迭代

采用高斯-赛德尔迭代法更新压力、速度场。

/*---------------------------高斯-赛德尔迭代-------------------------------
-------------------（采用松弛迭代,松弛因子为RXFACTOR)-----------------------
---1、P2 = P3 - RXFACTOR*(P3 - (-b/a*P3 - b/a*P2 - c/a*P3 - c/a*P2 - d/a))
---2、P3 = P2 - RXFACTOR*(P2 - (-b/a*P2 - b/a*P3 - c/a*P2 - c/a*P3 - d/a))
------------------------------------------------------------------------*/

#include "CavityFlow.h"

double GaussSeidel(double **U1, double **V1, double **P2, double **P3, int &itera1)
{
	int i, j, m;
	double a, b, c, d, e;
    double x, y;
	double k1, sum1;

    x = L/(BD-1), y = L/(BD-1);
	a = (T/(x*x)+T/(y*y))*2;
	b = -T/(x*x);
	c = -T/(y*y);

	for(m = 1; m <= SINGLESTEP; m++)
	{
		for(i = 2; i <= BD-1; i++)
		{
			for(j = 2; j <= BD-1; j++)
			{    
				if(i == 2 && j !=2 && j != BD-1)
			 	{   
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P3[i][j] - RXFACTOR * (P3[i][j] + b/(a+b)*P3[i+1][j] + c/(a+b)*P3[i][j+1] + c/(a+b)*P2[i][j-1] + d/(a+b));
				}
				else if(i == BD-1 && j! = 2 && j != BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P3[i][j] - RXFACTOR * (P3[i][j] + b/(a+b)*P2[i-1][j] + c/(a+b)*P3[i][j+1] + c/(a+b)*P2[i][j-1] + d/(a+b));
				}
				else if(j == 2 && i !=2 && i != BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P3[i][j] - RXFACTOR * (P3[i][j] + b/(a+c)*P3[i+1][j] + b/(a+c)*P2[i-1][j] + c/(a+c)*P3[i][j+1] + d/(a+c));
				}
				else if(j == BD-1 && i !=2 && i != BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P3[i][j] - RXFACTOR * (P3[i][j] + b/(a+c)*P3[i+1][j] + b/(a+c)*P2[i-1][j] + c/(a+c)*P2[i][j-1] + d/(a+c));
				}
				else if(i == 2 && j == 2)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P3[i][j] - RXFACTOR * (P3[i][j] + b/(a+b+c)*P3[i+1][j] + c/(a+b+c)*P3[i][j+1] + d/(a+b+c));
				}
				else if(i == 2 && j == BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P3[i][j] - RXFACTOR * (P3[i][j] + b/(a+b+c)*P3[i+1][j] + c/(a+b+c)*P2[i][j-1] + d/(a+b+c));
				}
				else if(i == BD-1 && j == 2)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P3[i][j] - RXFACTOR * (P3[i][j] + b/(a+b+c)*P2[i-1][j] + c/(a+b+c)*P3[i][j+1] + d/(a+b+c));
				}
				else if(i == BD-1 && j == BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P3[i][j] - RXFACTOR * (P3[i][j] + b/(a+b+c)*P2[i-1][j] + c/(a+b+c)*P2[i][j-1] + d/(a+b+c));
				}
				else
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P2[i][j] = P3[i][j] - RXFACTOR * (P3[i][j] + b/a*P3[i+1][j] + b/a*P2[i-1][j] + c/a*P3[i][j+1] + c/a*P2[i][j-1] + d/a);
				}
			}
		}

		for(i = 2; i <= BD-1; i++)
		{
			for(j = 2; j <= BD-1; j++)
			{    
				if(i == 2 && j !=2 && j != BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b)*P2[i+1][j] + c/(a+b)*P2[i][j+1] + c/(a+b)*P3[i][j-1] + d/(a+b));
				}
				else if(i == BD-1 && j !=2 && j != BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b)*P3[i-1][j] + c/(a+b)*P2[i][j+1] + c/(a+b)*P3[i][j-1] + d/(a+b));
				}
				else if(j == 2 && i !=2 && i != BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+c)*P2[i+1][j] + b/(a+c)*P3[i-1][j] + c/(a+c)*P2[i][j+1] + d/(a+c));
				}
				else if(j == BD-1 && i != 2 && i != BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+c)*P2[i+1][j] + b/(a+c)*P3[i-1][j] + c/(a+c)*P3[i][j-1] + d/(a+c));
				}
				else if(i == 2 && j == 2)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b+c)*P2[i+1][j] + c/(a+b+c)*P2[i][j+1] + d/(a+b+c));
				}
				else if(i == 2 && j == BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b+c)*P2[i+1][j] + c/(a+b+c)*P3[i][j-1] + d/(a+b+c));
				}
				else if(i == BD-1 && j == 2)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b+c)*P3[i-1][j] + c/(a+b+c)*P2[i][j+1] + d/(a+b+c));
				}
				else if(i == BD-1 && j == BD-1)
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/(a+b+c)*P3[i-1][j] + c/(a+b+c)*P3[i][j-1] + d/(a+b+c));
				}
				else
				{
                    d = 1/x*(U1[i][j]-U1[i-1][j])+1/y*(V1[i][j]-V1[i][j-1]);
					P3[i][j] = P2[i][j] - RXFACTOR * (P2[i][j] + b/a*P2[i+1][j] + b/a*P3[i-1][j] + c/a*P2[i][j+1] + c/a*P3[i][j-1] + d/a);
				}
			}
		}

		for(i = 2; i <= BD-1; i++)
		{
			for(j = 2; j <= BD-1; j++)
			{
				k1 = abs(P3[i][j]-P2[i][j]);
				sum1 = sum1+k1*k1;
			}
		}

		e = sqrt(sum1)/((BD-2)*(BD-2));
		if(e <= EPSILON2 && m != 1)
		{
			m = SINGLESTEP+1;
		}
		sum1 = 0;
		itera1++;
	}
	return e;
}

3.5 结果输出

主要用于计算结果输出，后处理软件采用Tecplot，所以保存文件类型为plt。

#include "CavityFlow.h"

void WriteResult(double **U, double **V, double **P, double Re)
{
    int i, j;
    int width = 0, height = 0;
    #ifdef GS
    string sstr = "result-"+to_string(int(Re))+"-"+to_string(BD)+"-"+to_string(T)+"-GS.plt";
    #else
    string sstr= "result-"+to_string(int(Re))+"-"+to_string(BD)+"-"+to_string(T)+"-JB.plt";
    #endif

    ofstream fout(sstr);

    fout<<"TITLE= \"Lid-Driven CavityFlow\""<<"\n";
    fout<<"VARIABLES=x,y,u,v,P"<<"\n";
    fout<<"ZONE I="<<BD<<",j="<<BD<<",F=POINT"<<"\n";

    for(j = 1; j <= BD; j++)
    {
        for(i = 1; i <= BD; i++)
        {
            fout<<width*(L/(BD-1))<<" "<<height*(L/(BD-1))<<" "<<(U[i][j]+U[i-1][j])/2<<" " <<(V[i][j]+V[i][j-1])/2<<" "<<P[i][j]<<"\n";
            width++;
            if(width == BD)
            {
                width = 0;
                height++;
            }
        }
    }
    fout<<endl;
    fout.close();
}

ofstream WriteResidual(double Re)
{
    #ifdef GS
	string istr= "tolerance-"+to_string(int(Re))+"-"+to_string(BD)+"-"+to_string(T)+"-GS.plt";
	#else
	string istr= "tolerance-"+to_string(int(Re))+"-"+to_string(BD)+"-"+to_string(T)+"-JB.plt";
	#endif

	ofstream fout(istr);
	fout<<"TITLE= \"Tolerance\""<<"\n";
    fout<<"VARIABLES=n,e,s"<<"\n";
    fout<<"ZONE t=line"<<"\n";

	return fout;
}

四、测试结果

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

Makefile如下：

SRC= LidDrivenCavityFlow.cpp Jacobi.cpp GaussSeidel.cpp OutputPlt.cpp
OBJ=$(SRC:.cpp=.o)
TARGET=LidDrivenCavityFlow

CC=g++
CFLAGS=-O2 -DGS

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


.PHONY:clean

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

编译：