格子Boltzamnn方法计算平板间方柱绕流流场速度压力分布

最新推荐文章于 2024-11-08 09:08:26 发布
最新推荐文章于 2024-11-08 09:08:26 发布
阅读量331 收藏 3
点赞数 1
文章标签： c++
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_45914943/article/details/123843718
版权 

/***********************************************
The LBE code is programed for isothermal incompressible flow using the D2Q9 with the bounce back boundary conditon.

for Square column flow around

**********************************************/
#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <fstream>
#include<fstream>
#include<sstream>
using namespace std;
const int N = 200;
const int M = 50;
const int DX = 20;
const int DY = 19;
const int H = 10;

const int N1 = (N + 1);
const int M1 = (M + 1);
const int jmax = M;
const int imax = N;
const int L = (M + 1);   // width of the cavity  
const int Q = 9;         // number of discrete velocities  
const double rho0 = 1.0; // initial density 
const double ux0 = 0.0;  // initial velocity component in x direction  
const double uy0 = 0.0;  // initial velocity component in y direction  
const double Re = 200.00;

const double Ma = 0.1;
const double RT = 1.0 / 3.0;
const double CS = sqrt(3 * RT);
const double Umax = CS * Ma;
//const double Umax = 0.1;
double UU;
int i, j, k;
int ex[Q] = { 0, 1, 0, -1, 0, 1, -1, -1, 1 };
int ey[Q] = { 0, 0, 1, 0, -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 f[M1][N1][Q], F[M1][N1][Q];
double ux[M1][N1], uy[M1][N1];
double u0[M1][N1], v0[M1][N1];                                      //error initial condition
double velocity[M1][N1], pressure[M1][N1], rho[M1][N1];             //calculate
double tau, err, err_out;                                           //Relaxtion time for BGK model
double C = 1.0;
void initial(void);
void evolution(void);
void boundary_P(void);
void calculate(void);
void dataoutput(void);
void contour(void);
double feq(double RHO, double U, double V, int k);
double error(void);
//================================================================
int main()
{
    int step, M2, N2;
    double err;
    M2 = M / 2;
    N2 = N / 2;

    step = 0;
    err = 1.0;
    tau = 3 * (double)H * Umax / Re + 0.5;


    clock_t start, finish;
    double time;
    start = clock();
    cout << "start calculate ..." << "\t  Re = " << Re << "\t tau = " << tau;
    cout << " \t  size: " << N << " * " << M << " \t    U = " << Umax << endl << endl;

    initial();

    while (err > 1.0e-6)
    {
        step++;
        evolution();
        calculate();
        boundary_P();
        if (step % 1000 == 0)
        {
            err = error();
            //printf("err = %e  \t ux_center = %e  \t uy_center = %e  \t step = %d\n", err, ux[M2][N2], uy[M2][N2], step);
            cout << scientific << "err = " << err << "\t  step = " << step << endl;
        }
        if (step >= 40000)
        {
            break;
        }
    }
    //dataoutput();
    contour();
    finish = clock();
    time = (float)(finish - start) / CLOCKS_PER_SEC;
    cout << "time = " << time << "s" << endl;
    return 0;
}
//================================================================

//1. 初始化
void initial()
{
    int i, j, k;
    for (j = 0; j <= jmax; j++)
    {
        for (i = 0; i <= imax; i++)
        {
            ux[j][i] = ux0;
            uy[j][i] = uy0;
            rho[j][i] = rho0;
            for (k = 0; k < Q; k++)
                f[j][i][k] = feq(rho[j][i], ux[j][i], uy[j][i], k);
        }
        double D = M;
        UU = (4 * Umax / D / D) * ((double)j * D - (double)j * j);
        ux[j][0] = UU;
    }
}

//2. 局部平衡分布函数
double feq(double RHO, double U, double V, int k)
{
    double euv, uv, re;
    euv = ex[k] * U + ey[k] * V;
    uv = U * U + V * V;
    re = w[k] * RHO * (1.0 + 3.0 * euv + 4.5 * euv * euv - 1.5 * uv);
    return re;
}

//3. 演化过程（碰撞迁移）
void evolution()
{
    int id, jd;
    for (j = 1; j < jmax; j++)//边界不参与碰撞
        for (i = 1; i < imax; i++)
        {
            if (i >= DX && i <= DX + H && j >= DY && j <= DY + H)
                continue;
            for (k = 0; k < Q; k++)
            {
                id = i - ex[k];
                jd = j - ey[k];
                F[j][i][k] = f[jd][id][k] + (feq(rho[jd][id], ux[jd][id], uy[jd][id], k) - f[jd][id][k]) / tau;
            }
        }
}

//4. 计算宏观量
void calculate()
{
    for (j = 1; j < jmax; j++)
    {
        for (i = 1; i < imax; i++)
        {
            if (i >= DX && i <= DX + H && j >= DY && j <= DY + H)
                continue;       //抛去方柱边界格点
            u0[j][i] = ux[j][i];//不需要收敛判断时可去除
            v0[j][i] = uy[j][i];
            rho[j][i] = 0.0;
            ux[j][i] = 0.0;
            uy[j][i] = 0.0;
            for (k = 0; k < Q; k++)
            {
                f[j][i][k] = F[j][i][k];
                rho[j][i] += f[j][i][k];    //求密度
                ux[j][i] += ex[k] * f[j][i][k];
                uy[j][i] += ey[k] * f[j][i][k];
            }
            ux[j][i] /= rho[j][i];
            uy[j][i] /= rho[j][i];
            velocity[j][i] = sqrt(ux[j][i] * ux[j][i] + uy[j][i] * uy[j][i]);   //计算宏观速度
            pressure[j][i] = rho[j][i] * C * C / 3;                             //计算压力
        }
    }

}

//5. 边界条件
void boundary_P()
{
    for (j = 0; j <= jmax; j++)//入口采用速度边界
    {
        i = 0;
        ux[j][i] = Umax;
        uy[j][i] = 0.0;
        rho[j][i] = (f[j][i][0] + f[j][i][2] + f[j][i][4] + 2.0 * (f[j][i][3] + f[j][i][6] + f[j][i][7])) / (1 - ux[j][i]);
        f[j][i][1] = f[j][i][3] + 2.0 * rho[j][i] * ux[j][i] / 3.0;
        f[j][i][5] = f[j][i][7] + rho[j][i] * ux[j][i] / 6.0 - 0.5 * (f[j][i][2] - f[j][i][4]) + rho[j][i] * uy[j][i] / 2.0;
        f[j][i][8] = f[j][i][6] + rho[j][i] * ux[j][i] / 6.0 + 0.5 * (f[j][i][2] - f[j][i][4]) - rho[j][i] * uy[j][i] / 2.0;
    }

    for (j = 0; j <= jmax; j++)//出口采用压力边界
    {
        i = N;
        ux[j][i] = 0;
        uy[j][i] = 0;
        rho[j][i] = 1.0;
        ux[j][i] = (f[j][i][0] + f[j][i][2] + f[j][i][4] + 2 * (f[j][i][1] + f[j][i][5] + f[j][i][8])) / rho[j][i] - 1.0;
        f[j][i][3] = f[j][i][1] - 2.0 * rho[j][i] * ux[j][i] / 3.0;
        f[j][i][7] = f[j][i][5] - rho[j][i] * ux[j][i] / 6.0 + 0.5 * (f[j][i][2] - f[j][i][4]);
        f[j][i][6] = f[j][i][8] - rho[j][i] * ux[j][i] / 6.0 - 0.5 * (f[j][i][2] - f[j][i][4]);
    }

    for (i = 1; i < N; i++)//上下边界采用标准反弹
    {
        j = 0;
        f[j][i][2] = f[j][i][4];
        f[j][i][5] = f[j][i][7];
        f[j][i][6] = f[j][i][8];
        j = M;
        f[j][i][4] = f[j][i][2];
        f[j][i][7] = f[j][i][5];
        f[j][i][8] = f[j][i][6];
    }

    for (i = DX + 1; i < DX + H; i++) //方柱上下边界采用标准反弹
    {
        j = DY + H;
        f[j][i][2] = f[j][i][4];
        f[j][i][5] = f[j][i][7];
        f[j][i][6] = f[j][i][8];
        j = DY;
        f[j][i][4] = f[j][i][2];
        f[j][i][7] = f[j][i][5];
        f[j][i][8] = f[j][i][6];
    }

    for (j = DY; j <= DY + H; j++) //方柱左右边界采用标准反弹
    {
        i = DX;
        f[j][i][3] = f[j][i][1];
        f[j][i][6] = f[j][i][8];
        f[j][i][7] = f[j][i][5];
        i = DX + H;
        f[j][i][1] = f[j][i][3];
        f[j][i][8] = f[j][i][6];
        f[j][i][5] = f[j][i][7];
    }
}

//6. 输出宏观量
void dataoutput()
{
    cout << "start output data ... " << endl;
    char n1[200];
    char n2[200];
    char n3[200];
    char n4[200];
    char n5[200];
    sprintf(n1, "x_%.0f.dat", Re);
    sprintf(n2, "y_%.0f.dat", Re);
    sprintf(n3, "ux_%.0f.dat", Re);
    sprintf(n4, "uy_%.0f.dat", Re);
    sprintf(n5, "rho_%.0f.dat", Re);

    int i, j;
    FILE* fp;

    fp = fopen(n1, "w+");
    for (i = 0; i <= N; i++) fprintf(fp, "%e \n", (i + 0.5) / L);
    fclose(fp);
    cout << "finish x.dat" << endl;

    fp = fopen(n2, "w+");
    for (j = 0; j <= M; j++) fprintf(fp, "%e \n", (j + 0.5) / L);
    fclose(fp);
    cout << "finish y.dat" << endl;

    fp = fopen(n3, "w");
    for (j = 0; j <= M; j++)
    {
        for (i = 0; i <= N; i++) fprintf(fp, "%e ", ux[j][i]);
        fprintf(fp, "\n");
    }
    fclose(fp);
    cout << "finish ux.dat" << endl;

    fp = fopen(n4, "w");
    for (j = 0; j <= M; j++)
    {
        for (i = 0; i <= N; i++) fprintf(fp, "%e ", uy[j][i]);
        fprintf(fp, "\n");
    }
    fclose(fp);
    cout << "finish uy.dat" << endl;

    fp = fopen(n5, "w");
    for (j = 0; j <= M; j++)
    {
        for (i = 0; i <= N; i++) fprintf(fp, "%e ", rho[j][i]);
        fprintf(fp, "\n");
    }
    fclose(fp);
    cout << "finish rho.dat" << endl;

}

//7. 输出云图
void contour()
{
    char c1[200];
    char c2[200];
    sprintf(c1, "Vector_%.0f .dat", Re);
    sprintf(c2, "Contour_%.0f.dat", Re);
    int i, j;
    //=======================================================
    ofstream fout;
    fout.open(c1, ios::app);
    fout << "TITLE = Vector   \nVariables = X Y U V  Pressure   Velocity \n" << "Zone " << "I = " << N1 << "  J = " << M1 << "   F = POINT" << endl;

    for (j = 0; j <= jmax; j++)
    {
        for (i = 0; i <= imax; i++)
        {
            fout << i << "\t" << j << "\t" << ux[j][i] << "\t" << uy[j][i] << "\t" << pressure[j][i] << "\t" << velocity[j][i] << endl;
        }
    }
    fout.close();
    cout << "Vector data finish" << endl;
    =======================================================
    //fout.open(c2, ios::app);
    //fout << "TITLE = Contour  \nVariables = X Y Velocity  \n" << "Zone " << "I = " << N1 << "  J = " << M1 << "   F = POINT" << endl;
    //for (j = 0; j <= jmax; j++)
    //{
    //    for (i = 0; i <= imax; i++)
    //    {
    //        fout << i << "\t" << j << "\t" << velocity[j][i] << endl;
    //    }
    //}
    //fout.close();
    //cout << "Contour data finish" << endl;
}

//8. 计算误差
double error()
{
    double err1, err2;
    err1 = err2 = 0.0;
    for (j = 1; j < jmax; j++)
        for (i = 1; i < imax; i++)
        {
            if (i >= DX && i <= DX + H && j >= DY && j <= DY + H)
                continue;
            err1 += sqrt((ux[j][i] - u0[j][i]) * (ux[j][i] - u0[j][i]) + (uy[j][i] - v0[j][i]) * (uy[j][i] - v0[j][i]));
            err2 += sqrt(ux[j][i] * ux[j][i] + uy[j][i] * uy[j][i]);
        }
    return err1 / (err2 + 1e-30);
}

【论文阅读】 AMNN: Attention-Based Multimodal Neural Network Model for Hashtag Recommendation
qq_43146737的博客
11-15 3173
这里写自定义目录标题摘要引言相关的工作我们的模型 摘要 在真实的社交网络中，标签被广泛应用于理解单个微博的内容。然而，用户在发布微博时并不总是主动添加话题标签，因此需要投入大量精力进行自动的话题标签推荐。作为一种新趋势，用户不再仅仅发布文本，而是更愿意分享多模态数据，如图像。针对这种情况，我们提出了一种基于注意力的多模态神经网络模型(AMNN)来学习多模态微博的表示并推荐相关的标签。在本文中，我们将hashtag推荐任务转换为一个序列生成问题。然后，我们提出了一种混合神经网络方法来提取文本和图像的特征，并将
Python-LBM(格子玻尔兹曼)程序源码实例分析—圆柱绕流篇
热门推荐
weixin_58913471的博客
06-12 1万+
初次学习LBM计算方法，找到一个比较优秀的用python语言编写的圆柱绕流的实例，对每段代码详细添加了注释，帮助自己总结，也为初学者提供一点帮助。 先放一张结果图像：
LBM方柱绕流C++.docx
10-29
采用格子Boltzmann方法（LBM）对流体力学经典问题方柱绕流进行分析计算，得出的数据可用后处理软件进行处理。本代码仅供学习之用，对于新手较为友好，如有问题，请指正。欢迎大家在评论区发言
有限元法计算二维圆柱绕流问题——Python代码实现
weixin_52046577的博客
05-11 1687
选取流函数Ψ为变量，对拉普拉斯方程进行求解（右边界为自然边界条件，其余边界为本质边界条件）；网格数据文件的生成暂时不在本文中详述。
平板间方柱绕流的格子Boltzmann方法模拟 (2010年)
05-17
为研究平板间方柱绕流上下平板对置于其中的方柱绕流所产生的影响，采用格子Boltzmann方法对二维平板间低雷诺数(Re=100)方柱绕流问题进行了数值模拟研究．分析了3种不同阻塞比下，平板边壁对方柱的升、阻力系数、...
如何利用C++实现格子Boltzmann方法（LBM）来模拟方柱绕流问题，并使用后处理软件进行流场分析？
10-30
为了解决如何利用C++实现格子Boltzmann方法（LBM）来模拟方柱绕流问题并进行流场分析的问题，可以参考以下步骤进行操作：首先，确保你已经安装了用于编写C++的开发环境，例如Visual Studio或者GCC。接着，参考《LBM...
如何使用格子Boltzmann方法（LBM）在C++中编写代码来模拟方柱绕流问题？
10-30
在流体力学领域，格子Boltzmann方法（LBM）是一种有效的数值模拟手段，尤其适用于复杂流场的计算。为了帮助你掌握如何在C++中使用LBM进行方柱绕流问题的模拟，我建议你参考这本资料《LBM方法求解方柱绕流问题的C++...
能否指导我如何利用C++和格子Boltzmann方法（LBM）编写代码来模拟方柱绕流问题，并指出哪些后处理软件适合分析模拟结果？
10-30
格子Boltzmann方法（LBM）是计算流体力学中的一种有效技术，特别适用于复杂的流动问题，如方柱绕流。在C++中编写LBM代码来模拟方柱绕流涉及以下步骤： 参考资源链接：[LBM方法求解方柱绕流问题的C++代码实现]...
c语言求解非线性方程组_篇外6、无粘平板扰流的求解
weixin_39764212的博客
12-25 173
首先声明，本篇内容改编自本人本科三年级计算流体力学课程设计。由于现在的主要精力在准备考研，所以可能需要几天才会想起来还有这么一个专栏，为了能尽快更（shui）篇目，就决定分别将自己大三下半学期的计算流体力学课程设计和大四上学期（也就是当前学期）的流体力学课程设计的内容放上来。内容整体上比较初级，不是很难的内容，但是为了能让最初的初学者都能看懂所以比较碎叨（因为我当初自学时就因为看不懂受到了很大挫折...
简单代码=基于LBM二维圆柱绕流代码
04-07
二维圆柱绕流简单LBM代码
网格问题代码
weixin_41655489的博客
04-11 860
#include &lt;iostream&gt;                                      //声明头文件#include &lt;cmath&gt;#include &lt;cstdlib&gt;#include &lt;iomanip&gt;#include &lt;fstream&gt;#include &lt;sstream&gt;#include &lt...
格子玻尔兹曼法学习记录（附MATLAB画图源程序）
weixin_37783345的博客
03-30 9140
感谢群里的朋友们提供帮助！还是老样子，有啥问题Feel free to tell us~毕竟群众力量大嘛~格子玻尔兹曼救星QQ群：293267908。 流体计算领域中，LBM还是个比较新的思想，最近宝宝正在尝试性的进行研究。首先说一下教材吧。我在某宝买了中文版，杨大勇翻译的，中国工信出版。 还可以看英文原版哈，英文教材链接：https://pan.baidu.c...
LBM盖顶驱动流C++代码
diloudiao5379的博客
03-05 2143
#include<iostream> #include<cmath> #include<cstdlib> #include<iomanip> #include<fstream> #include<sstream> #include<string> using namespace st...
格子玻尔兹曼机(Lattice Boltzmann Method)系列4：LBM实例之方腔驱动流
weixin_43796562的博客
02-23 3977
问题描述 在一个正方形的空腔中，顶盖以U的速度被突然启动，导致产生了一个在边角处会有小的漩涡出现的流场。 左、右、下边界使用了回弹性边界条件，而上边界使用了上一篇文章中的Von Neumann速度边界条件。 下面贴出来具体的代码，环境是：VS2013 + OpenMP +OpenGL。如果使用VS的话，不用另外下载任何东西，只需要简单的几步设置就可以配置出OpenMP，具体打开开关的步骤直接百度...
编译器、IDE对C/C++新标准的支持
最新发布
编程语言技巧经验分享
11-08 747
很遗憾的是，微软VC++对C99的支持直到VS2015才基本完全支持，在VS2012之前，仅仅支持C89. 也正是2015年之前，这么多人骂微软并纷纷转移到其它C语言编译器的原因。VS2019开始才正式宣布支持C11/C17.对于C++标准，VS2013开始支持C++11，VS2017支持C++17, VS2019支持C++20, 显然，微软更偏爱C++.
C++文件的读取操作
m0_67677309的博客
11-04 157
fstream 对文件的读写操作，继承了ofstream\ifstream类的功能。1. ofstream 对文件的写操作，继承了ostream类的功能。ifstream 对文件的读操作，继承了istream类的功能。对文件操作，打开才能写。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Yvie_Bae

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值