应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

Znn3bq.jpeg
小木虫论坛-学术科研互动平台 » 专业学科区 » 机械 » 流体机械 » 求助格子boltzmann方法模拟平板间气体流动
61/1返回列表
查看: 1128  |  回复: 5
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

花开时节0931

铁虫 (初入文坛)

[求助] 求助格子boltzmann方法模拟平板间气体流动 已有1人参与

空气以声速从平板间进入，平板厚度10微米，想采用格子boltzmann模拟气体流动，没学过C语言，边学边改。问题还是很大，毕业论文着急，求指导。不胜感激。
#include <cmath>
#include <iostream>
#include <cstdlib>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
using namespace std;
const int Q=9; //D2Q9模型
const int NX=120;//X方向
const int NY=40;//Y方向
//const double U=0.1; //顶盖速度
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[NX+1][NY+1],u[NX+1][NY+1][2],u0[NX+1][NY+1][2],f[NX+1][NY+1][Q],F[NX+1][NY+1][Q];
int i,j,k,ip,jp,n;
double c,Re,dx,dy,Lx,Ly,dt,rho0,P0,tau_f,niu,error,Kn;
void init();
double feq(int k,double rho,double u[2]);
double rhosum[NX][NY] = {0};
double Psum[NX][NY] = {0};

void evolution();
void output(int m) ;
void Error();
int main()
{
    using namespace std;
    init();
    for(n=0; ;n++)
    {
        evolution();
        if (n%1000==0)
        {
             Error();
             cout<<"the "<<n<<" th computation result:"<<endl<<
             "The u,v of point (NX/2,NY/2)is:" <<setprecision(6)
             <<u[NX/2][NY/2][0]<<","<<u[NX/2][NY/2][1]<<endl;
             cout<<"The max relative error of uv is: "
             <<setiosflags(ios::scientific)<< error <<endl;
             if(n>=100)
             {
                if(n%1000==0) output(n);
                if(error<1.0e-12) break;
             }
        }
    }
    return 0;
}
void init()
{
  dx=1.0;
  dy=1.0;
  P0=0.5;
  Lx=dx*double(NY);
  Ly=dy*double(NX);
  dt=dx;
  c=dx/dt;//1.0
  rho0=1.03;
  Kn=0.5;
  tau_f=8;
//tau_f=P0*Lx*Kn*sqrt(6/3.14)+0.5;
  //Re=8;
// niu=0.1*Lx/Re;
//tau_f=3.0*niu+0.5;

  std::cout<<"tau_f="<<tau_f<<endl;
      
  for(i=0;i<=NX;i++)//初始化
  for(j=0;j<=NY;j++)
  {
        u[j][0]=0;
        u[j][1]=0;
        rho[j]=rho0;
        u[0][j][0]=0.57;
        for (k=0;k<Q;k++)
        {
            f[j][k]=feq(k,rho[j],u[j]);

        }
  }
}
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;
}

void evolution()
{
for(i=1;i<NX;i++)//演化
for(j=1;j<NY;j++)
for(k=0;k<Q;k++)
{
      ip=i-e[k][0];
      jp=j-e[k][1];
      F[j][k]=f[ip][jp][k]+(feq(k,rho[ip][jp],
      u[ip][jp])-f[ip][jp][k])/tau_f;
}
         
for(i=1;i<NX;i++)//计算宏观量
  for(j=1;j<NY;j++)
  {
       u0[j][0]=u[j][0];
       u0[j][1]=u[j][1];
       rho[j]=0;
       u[j][0]=0;
       u[j][1]=0;
       for(k=0;k<Q;k++)
       {
           f[j][k]=F[j][k];
           rho[j]+=f[j][k];
           u[j][0]+=e[k][0]*f[j][k];
           u[j][1]+=e[k][1]*f[j][k];
                   rhosum[j] += f[j][k];                  
       }
           Psum[j] = rhosum[j] * c * c / 3;
       u[j][0]/=rho[j];
       u[j][1]/=rho[j];
   }
           
//边界处理
//左右边界  
  

for(i=0;i<=NX;i++)//上下边界
  for(k=0;k<Q;k++)
  {
       rho[0]=rho[1];
       f[0][k]=feq(k,rho[0],u[0])+f[1][k]
       -feq(k,rho[1],u[1]);
       rho[NY]=rho[NY-1];
       u[NY][0]=0;
       u[NY][1]=0;
           u[0][0]=0;
       u[0][1]=0;

       f[NY][k]=feq(k,rho[NY],u[NY])+f[NY-1]
       [k]-feq(k,rho[NY-1],u[NY-1]);
  }
}

void output(int m)//输出
{
      ostringstream name;
      name<<"cavity_"<< m<<".dat";
      ofstream out(name.str().c_str());
      out<<"Title=\"LBM Lid Driven Flow\"\n"<<
      "VARIABLES=\"X\",\"Y\", \"U\",\"V\",\"rhosum\",\"Psum\"\n"
      <<"ZONE T=\"BOX\",I="<<NX+1<<",J="<<NY+1<<",F=POINT"<<endl;
      for(j=0;j<=NY;j++)
      for(i=0;i<=NX;i++)
      {
        out<<
        double(i)/Lx<<","<<double(j)/Ly<<","<<u[j][0]<<","<<u[j][1]<<","<<rhosum[j]<<","<<Psum[j]<<"\n"
        <<endl;
      }
}
void Error()
{
     double temp1,temp2;
     temp1=0;
     temp2=0;
     for(i=1;i<NX;i++)
        for(j=1;j<NY;j++)
        {
             temp1+=(
             (u[j][0]-u0[j][0])*(u[j][0]-u0[j][0])
             +(u[j][1]-u0[j][1])*(u[j][1]-u0[j][1]));
              temp2+=
             (u[j][0]*u[j][0]+u[j][1]*u[j][1]);
                        
        }
    temp1=sqrt(temp1);
    temp2=sqrt(temp2);
    error=temp1/(temp2);
}
回复此楼

» 猜你喜欢
1楼 2015-11-22 22:29:34
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老大读书

铜虫 (正式写手)
有什么问题

发自小木虫Android客户端
回复此楼
2楼2015-11-24 23:28:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

独孤694

至尊木虫 (知名作家)

【答案】应助回帖

感谢参与，应助指数 +1
你是根据何雅玲书里的程序改的？没细看，首先她程序里是不可压模型，速度要保证小于0.3马赫，而你的问题是可压流，目前最常用的可压模型是多速模型吧

发自小木虫Android客户端
赞一下 回复此楼
人就得逼自己一把
3楼2015-11-24 23:39:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

孤独的漫步

木虫 (小有名气)
为啥里面程序的‘ ’都没了。。。
赞一下 回复此楼
你必须非常努力，才能看起来毫不费力
4楼2015-11-29 19:43:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

花开时节0931

铁虫 (初入文坛)
引用回帖:
2楼: Originally posted by 老大读书 at 2015-11-24 23:28:52
有什么问题

想请问下有没有压力驱动的程序参考，现在只会速度驱动
赞一下 回复此楼
5楼2015-12-29 10:43:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

花开时节0931

铁虫 (初入文坛)
引用回帖:
3楼: Originally posted by 独孤694 at 2015-11-24 23:39:03
你是根据何雅玲书里的程序改的？没细看，首先她程序里是不可压模型，速度要保证小于0.3马赫，而你的问题是可压流，目前最常用的可压模型是多速模型吧

是的 小白 想做压力驱动 但是只知道边界上改 初始设置不知道怎么设置 请问有没有想过压力驱动程序参考下
赞一下 回复此楼
6楼2015-12-29 10:45:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 花开时节0931 的主题更新
61/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 293求调剂 +10 我爱高数高数爱� 2026-04-12 10/500 2026-04-12 11:15 by 一切OK
[考研] 电气专硕320求调剂 +6 小麻子111 2026-04-10 6/300 2026-04-12 10:54 by lemon6009
[考研] 0854调剂 +8 长弓傲 2026-04-11 9/450 2026-04-11 21:44 by 蓝云思雨
[考研] 295求调剂 +6 ?要上岸? 2026-04-05 7/350 2026-04-11 19:02 by laoshidan
[考研] 化工求调剂！ +35 RichLi_ 2026-04-06 35/1750 2026-04-11 11:02 by zhq0425
[考研] 生物学调剂 +8 小冉要努力 2026-04-10 9/450 2026-04-11 10:22 by wwj2530616
[考研] 275求调剂 +9 1624447980 2026-04-08 10/500 2026-04-11 10:20 by Delta2012
[考研] 0854调剂 +4 长弓傲 2026-04-09 4/200 2026-04-11 09:18 by 猪会飞
[考研] 263能源动力专硕求调剂 +4 加大号饭盒袋 2026-04-10 4/200 2026-04-10 20:52 by gong120082
[考研] 一志愿京区985，085401电子信息，本科电子信息 +3 阳光开朗的男孩 2026-04-10 3/150 2026-04-10 16:29 by sophia_93
[考研] 301求调剂 +5 149. 2026-04-10 5/250 2026-04-10 15:45 by 柴小白
[考研] 336材料与化工085600求调剂 +21 水星记infp 2026-04-05 24/1200 2026-04-10 15:28 by luoyongfeng
[考研] 一志愿沪9，326生物学求相关专业调剂 +4 刘墨墨 2026-04-09 4/200 2026-04-10 12:07 by pengliang8036
[考研] 材料化工总分334求调剂 +16 Riot2025 2026-04-08 17/850 2026-04-09 20:19 by maddjdld
[论文投稿] 求助文献原文 10+3 18500821399 2026-04-08 3/150 2026-04-09 16:56 by 北京莱茵润色
[考研] 化学工程与技术专业一志愿哈工程 291分B区 国家级大创负责人 有一作论文 +13 Emmy~ 2026-04-09 13/650 2026-04-09 14:47 by only周
[考研] 生物学308分求调剂（一志愿华东师大） +13 相信必会光芒万� 2026-04-06 16/800 2026-04-09 13:54 by 徐良白眉大侠
[考研] 270求调剂 +3 031127 2026-04-06 4/200 2026-04-08 21:00 by 逆水乘风
[考研] 071000生物学，一志愿深圳大学296分，求调剂 +12 TIckLw 2026-04-06 13/650 2026-04-07 20:34 by lijunpoly
[考研] 362求调剂一志愿中国石油大学 +4 我要考大 2026-04-06 6/300 2026-04-06 14:11 by 无际的草原
信息提示
关闭
请填处理意见
关闭