目录
一、物理模型
顶盖驱动方腔流(Lid-driven Cavity Flow)是经典的计算流体力学(CFD)问题,也是验证流体力学问题计算方法的基准案例之一。二维顶盖驱动方腔流问题的物理模型由一个二维方形空腔构成,其上壁面具有水平方向的速度,其它三个壁面为静止壁面(速度为0),模型如图1所示。流体力学实验表明,该形式的运动会产生一种流动现象,特征是:空腔中心会形成一个大涡,空腔角落形成一些小涡,而雷诺数的大小将会影响涡的大小和数量。流动特征如图2所示。
二、数学模型
2.1 流动控制方程
采用不可压缩定常流Navier-Stokes方程:
将公式(1)中速度、压力、长度参数无量纲化后得到不可压缩定常流Navier-Stokes方程:
将公式(2)写成x和y方向的二维方程组,可得:
2.2 离散方法
1、求解域离散
计算域为方形,边长,采用二维均匀网格划分,将计算域离散为
的均匀网格,空间步长为
,离散后的网格模型如下图所示(n=200)。

2、控制方程离散
模型的边界条件为:
使用交错网格,将U和V的边界条件视为虚拟节点与内部一层节点之和的1/2等于0,即虚拟节点与内部节点互为相反数。P的边界条件在网格中表现为第一层和第二层(或最后一层和倒数第二层)相等。
差分方程:
将假设的U和V代入差分方程(4a)和(4b),得到当前的U*和V*,如下式:
其中:
2.3 迭代方法
迭代方法分别采用雅各比(Jacobi)迭代和高斯-赛德尔(Gauss-Seidel)迭代来求解,这两种迭代方法的具体理论推导请参考博主的另一篇博文:↓传送门↓
偏微分方程算法之迭代法(番外篇)_迭代法解偏微分方程-优快云博客https://blog.youkuaiyun.com/L_peanut/article/details/138302224 将U*和V*代入公式(6),可得到P',利用雅各比迭代或高斯-赛德尔迭代求解方程组:
其中:
将P'代入中,得到当前时刻准确的P值,再将P'代入公式(7a)和(7b),得到当前时刻准确的U和V值。
再将当前时刻准确的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)
编译:
4.1 顶盖速度U=0.05,Re=500
计算结果:
X方向速度分布:
流线图:
4.2 顶盖速度U=0.3,Re=3000
计算结果:
X方向速度分布:
流线图:
4.3 顶盖速度U=1.0,Re=10000
计算结果:
X方向速度分布:
流线图:
由三种不同雷诺数计算结果可知,当Re较小时,方腔内部未形成局部小涡;当Re逐渐增大时,方腔内部形成局部小涡。当Re为5000时,方腔底部两侧均形成一个小涡;Re达到10000时,涡的数量增加,位置发生变化,说明流体湍动程度增大。