//Program By Yuanqing XU-2022-11-10
#include “stdafx.h”
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
using namespace std;
//基本参数，需要从软件前端设置
const int Lx = 200; //x方向网格数 （参数1）
const int Ly = 50; //y方向网格数 （参数2）
const int MaxT = 3000; //迭代次数 （参数3）
const int OutN = 200; //输出间隔 （参数4）
const double Uin = 0.1; //入口速度 （参数5）
double Re = 30; // 雷诺数 (参数6)
const char *const PathD = “..\\DataOut\\U_”; // 数据存放路径 (参数7)
const double pi = 3.14159265;
//LBM参数
const int Q = 9; //D2Q9模型
double rho0 = 1.0;
double nu = Uin*Ly / Re;
double tau = 3.0*nu + 0.5;
//变量初始化
int i, j, k, ip, jp, tStep, ii, jj;
int e[Q][2] = { { 0, 0 }, { 1, 0 }, { 0, 1 }, { -1, 0 }, { 0, -1 }, { 1, 1 }, { -1, 1 }, { -1, -1 }, { 1, -1 } };
double w[Q] = { 4.0 / 9, 1.0 / 9, 1.0 / 9, 1.0 / 9, 1.0 / 9, 1.0 / 36, 1.0 / 36, 1.0 / 36, 1.0 / 36 };
double rho[Ly + 1][Lx + 1], F[Ly + 1][Lx + 1][Q], FEQM[Ly + 1][Lx + 1][Q], U[Ly + 1][Lx + 1][2];
double F0[Ly + 1][Lx + 1][Q];
//函数初始化
double feq(int k, double rho, double u[2]);
void Initialization();
void Outputdata(int m);
//主程序
int main()
{
using namespace std;
Initialization();
for (tStep = 1; tStep <= MaxT; tStep++)
{
//1- 迁移（不包含边界）
for (i = 1; i < Lx; i++)
for (j = 1; j < Ly; j++)
for (k = 0; k < Q; k++)
{
{
{
ii = i – e[k][0]; jj = j – e[k][1];
F[j][i][k] = F0[jj][ii][k];
}
}
}
//2- 计算宏观量
for (i = 1; i < Lx; i++)
for (j = 1; j < Ly; j++)
{
{
rho[j][i] = 0;
U[j][i][0] = 0;
U[j][i][1] = 0;
for (k = 0; k < Q; k++)
{
rho[j][i] += F[j][i][k];
U[j][i][0] += e[k][0] * F[j][i][k];
U[j][i][1] += e[k][1] * F[j][i][k];
}
U[j][i][0] /= rho[j][i];
U[j][i][1] /= rho[j][i];
}
}
//左右宏观边界
for (j = 1; j < Ly; j++)
{
U[j][0][0] = Uin; U[j][0][1] = 0;
U[j][Lx][0] = U[j][Lx – 1][0]; U[j][Lx][1] = 0;
rho[j][0] = rho[j][1]; rho[j][Lx] = rho0;
}
//上下壁面
for (i = 0; i <= Lx; i++)
{
U[0][i][0] = 0; U[0][i][1] = 0;
U[Ly][i][0] = 0; U[Ly][i][1] = 0;
rho[0][i] = rho[1][i];
rho[Ly][i] = rho[Ly – 1][i];
}
//4- 碰撞（不包含边界）
for (i = 1; i < Lx; i++)
for (j = 1; j < Ly; j++)
for (k = 0; k < Q; k++)
{
{
{
F0[j][i][k] = F[j][i][k] – 1 / tau*(F[j][i][k] – feq(k, rho[j][i], U[j][i]));
}
}
}
//5- 微观边界处理（不需要判断类型和方向，处理都一样）
//左右微观边界
for (j = 1; j < Ly; j++)
for (k = 0; k < Q; k++)
{
{
F0[j][0][k] = feq(k, rho[j][0], U[j][0]) + F0[j][1][k] – feq(k, rho[j][1], U[j][1]); //入口
F0[j][Lx][k] = feq(k, rho[j][Lx], U[j][Lx]) + F0[j][Lx – 1][k] – feq(k, rho[j][Lx – 1], U[j][Lx – 1]); //出口
}
}
//上下微观边界
for (i = 1; i < Lx; i++)
for (k = 0; k < Q; k++)
{
{
F0[0][i][k] = feq(k, rho[0][i], U[0][i]) + F0[1][i][k] – feq(k, rho[1][i], U[1][i]); //下边界
F0[Ly][i][k] = feq(k, rho[Ly][i], U[Ly][i]) + F0[Ly – 1][i][k] – feq(k, rho[Ly – 1][i], U[Ly – 1][i]); //上边界
}
}
//6 输出数据和显示
if (tStep % OutN == 0)
{
cout << “入口中点速度：” << U[25][1][0] << ” NO.” <<tStep<< “—Steps” << endl;
Outputdata(tStep);
}
}
return 0;
}
// 主程序结束~~~~~~~~~~~
//子程序与函数======================================================================
//子程序1 初始化函数
void Initialization()
{
for (i = 0; i <= Lx; i++)
for (j = 0; j <= Ly; j++)
{
U[j][i][0] = 0; U[j][i][1] = 0;
rho[j][i] = 1; U[j][0][0] = Uin;
for (k = 0; k < Q; k++)
{
F[j][i][k] = feq(k, rho[j][i], U[j][i]);
F0[j][i][k] = feq(k, rho[j][i], U[j][i]);
}
}
U[0][0][0] = 0;
U[Ly][0][0] = 0;
}
////子程序2 计算输出速度场函数
void Outputdata(int m)
{
ostringstream name;
name << PathD << m << “.dat”; //数据存储地址
ofstream out(name.str().c_str());
out << “Title=\”LBM Channel Flow\”\n”
<< “VARIABLES=\”X\”,\”Y\”,\”U\”,\”V\”\n”
<< “ZONE T= \”BOX\”, I= “
<< Lx + 1 << “, J=” << Ly + 1 << “, F=POINT” << endl;
for (j = 0; j <= Ly; j++)
for (i = 0; i <= Lx; i++)
{
out << double(i) << ” ” << double(j) << ” “
<< U[j][i][0] << ” ” << U[j][i][1] << endl;
}
}
//函数1 计算平衡态分布函数
double feq(int k, double rho, double u[2])
{
double eu, uv, feq;
eu = (e[k][0] * u[0] + e[k][1] * u[1]);
uv = (u[0] * u[0] + u[1] * u[1]);
feq = w[k] * rho*(1.0 + 3.0*eu + 4.5*eu*eu – 1.5*uv);
return feq;
}