花開時(shí)節(jié)0931

鐵蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 21
帖子: 19
在線: 5.8小時(shí)
蟲號: 3159261
注冊: 2014-04-24

[求助] 求助格子boltzmann方法模擬平板間氣體流動已有1人參與

空氣以聲速從平板間進(jìn)入，平板厚度10微米，想采用格子boltzmann模擬氣體流動，沒學(xué)過C語言，邊學(xué)邊改。問題還是很大，畢業(yè)論文著急，求指導(dǎo)。不勝感激。
#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])//計(jì)算平衡態(tài)分布函數(shù)
{
   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++)//計(jì)算宏觀量
  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);
}

回復(fù)此樓

» 猜你喜歡

求求2Cr13的激光沖擊有限元模擬J-C本構(gòu)參數(shù) 已經(jīng)有1人回復(fù)
上海交大材料學(xué)院招26年秋“申請-考核制”生物材料方向普博生1名 (截至4月10日) 已經(jīng)有8人回復(fù)
金屬材料論文潤色/翻譯怎么收費(fèi)? 已經(jīng)有254人回復(fù)
上海交大材料學(xué)院招26年秋“申請-考核制”生物材料方向普博生1名 (截至4月10日) 已經(jīng)有7人回復(fù)
上海交大材料學(xué)院招26年秋“申請-考核制”生物材料方向普博生1名 (截至4月10日) 已經(jīng)有6人回復(fù)
中國科學(xué)院-材料計(jì)算小組尋求2026級研究生（包括調(diào)劑生/聯(lián)合培養(yǎng)等）已經(jīng)有132人回復(fù)
天津高校航空材料先進(jìn)制造技術(shù)與裝備課題組擬接收調(diào)劑碩士研究生5-6名。已經(jīng)有0人回復(fù)
0858專碩.339分求調(diào)劑已經(jīng)有0人回復(fù)
294 英二數(shù)二物化求調(diào)劑已經(jīng)有6人回復(fù)
找金屬材料類期刊已經(jīng)有2人回復(fù)
西華大學(xué)材料學(xué)院表面科學(xué)與工程技術(shù)科研團(tuán)隊(duì)2026年招收研究生已經(jīng)有0人回復(fù)

1樓 2015-11-22 22:29:34

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

老大讀書

銅蟲 (正式寫手)

應(yīng)助: 5 (幼兒園)
金幣: 165.4
散金: 305
紅花: 2
帖子: 560
在線: 338.5小時(shí)
蟲號: 2471699
注冊: 2013-05-19
專業(yè): 電化學(xué)分析

有什么問題

發(fā)自小木蟲Android客戶端

回復(fù)此樓

2樓2015-11-24 23:28:52

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

獨(dú)孤694

至尊木蟲 (知名作家)

應(yīng)助: 231 (大學(xué)生)
金幣: 6733.6
散金: 9022
紅花: 41
沙發(fā): 23
帖子: 5601
在線: 1898.6小時(shí)
蟲號: 2134245
注冊: 2012-11-18
性別: GG
專業(yè): 多相流熱物理學(xué)

【答案】應(yīng)助回帖

感謝參與，應(yīng)助指數(shù) +1

你是根據(jù)何雅玲書里的程序改的？沒細(xì)看，首先她程序里是不可壓模型，速度要保證小于0.3馬赫，而你的問題是可壓流，目前最常用的可壓模型是多速模型吧

發(fā)自小木蟲Android客戶端

贊一下

回復(fù)此樓

人就得逼自己一把

3樓2015-11-24 23:39:03

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

孤獨(dú)的漫步

木蟲 (小有名氣)

應(yīng)助: 6 (幼兒園)
金幣: 4318.8
散金: 677
紅花: 2
帖子: 87
在線: 352小時(shí)
蟲號: 2591375
注冊: 2013-08-11
性別: GG
專業(yè): 傳熱傳質(zhì)學(xué)

為啥里面程序的‘ ’都沒了。。。

贊一下

回復(fù)此樓

你必須非常努力，才能看起來毫不費(fèi)力

4樓2015-11-29 19:43:42

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

花開時(shí)節(jié)0931

鐵蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 21
帖子: 19
在線: 5.8小時(shí)
蟲號: 3159261
注冊: 2014-04-24

引用回帖:

2樓: Originally posted by 老大讀書 at 2015-11-24 23:28:52
有什么問題

想請問下有沒有壓力驅(qū)動的程序參考，現(xiàn)在只會速度驅(qū)動

贊一下

回復(fù)此樓

5樓2015-12-29 10:43:57

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

花開時(shí)節(jié)0931

鐵蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 21
帖子: 19
在線: 5.8小時(shí)
蟲號: 3159261
注冊: 2014-04-24

引用回帖:

3樓: Originally posted by 獨(dú)孤694 at 2015-11-24 23:39:03
你是根據(jù)何雅玲書里的程序改的？沒細(xì)看，首先她程序里是不可壓模型，速度要保證小于0.3馬赫，而你的問題是可壓流，目前最常用的可壓模型是多速模型吧

是的小白想做壓力驅(qū)動但是只知道邊界上改初始設(shè)置不知道怎么設(shè)置請問有沒有想過壓力驅(qū)動程序參考下

贊一下

回復(fù)此樓

6樓2015-12-29 10:45:44

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

最具人氣熱帖推薦 [查看全部]	作者	回/看	最后發(fā)表

[考研] 一志愿太原理工安全工程300分，求調(diào)劑 +5	0857求調(diào)劑. 2026-03-24	6/300	2026-03-28 22:04 by zhq0425
[考研] 283求調(diào)劑 +3	A child 2026-03-28	3/150	2026-03-28 15:41 by ms629
[考研] 312，生物學(xué)求調(diào)劑 +3	小譯同學(xué)abc 2026-03-28	3/150	2026-03-28 15:32 by 落�？伤�
[考研] 085600，材料與化工321分求調(diào)劑 +9	大饞小子 2026-03-28	9/450	2026-03-28 14:56 by 神馬都不懂
[考研] 346求調(diào)劑一志愿070303有機(jī)化學(xué) +3	蘿卜燉青菜 2026-03-28	3/150	2026-03-28 14:11 by 唐沐兒
[考研] 266分，求材料冶金能源化工等調(diào)劑 +7	哇呼哼呼哼 2026-03-27	9/450	2026-03-28 12:22 by zllcz
[考研] 286求調(diào)劑 +4	丟掉懶惰 2026-03-27	7/350	2026-03-28 08:07 by baoball
[考研] 086502化學(xué)工程342求調(diào)劑 +6	阿姨復(fù)古不過 2026-03-27	6/300	2026-03-28 07:06 by wangy0907
[考研] 一志愿南師大0703化學(xué) 275求調(diào)劑 +4	Ripcord上岸 2026-03-27	4/200	2026-03-27 17:00 by zhyzzh
[考研] 279 分求調(diào)劑 +4	睡個(gè)好覺_16 2026-03-24	4/200	2026-03-27 15:05 by 醉在風(fēng)里
[考研] 284求調(diào)劑 +11	junqihahaha 2026-03-26	12/600	2026-03-27 04:37 by wxiongid
[考研] 342求調(diào)劑 +3	加油a李zs 2026-03-26	3/150	2026-03-27 00:29 by wxiongid
[考研] 求調(diào)劑 +8	Auroracx 2026-03-22	8/400	2026-03-26 19:55 by 不吃魚的貓
[考研] 一志愿河工大 081700 276求調(diào)劑 +4	地球繞著太陽轉(zhuǎn) 2026-03-23	4/200	2026-03-26 14:27 by zzll406
[考研] 化學(xué)調(diào)劑一志愿上海交通大學(xué)336分-本科上海211 +4	小魚愛有機(jī) 2026-03-25	4/200	2026-03-26 10:19 by aa331100
[考研] 332求調(diào)劑 +6	032500 2026-03-25	6/300	2026-03-25 22:45 by 418490947
[考研] 各位老師您好：本人初試372分 +5	jj涌77 2026-03-25	6/300	2026-03-25 14:15 by mapenggao
[考研] 282求調(diào)劑 +3	wcq131415 2026-03-24	3/150	2026-03-25 12:16 by userper
[考研] 生物學(xué)學(xué)碩求調(diào)劑 +7	小羊睡著了? 2026-03-23	10/500	2026-03-25 02:24 by 清風(fēng)拂揚(yáng)。 m
[考研] 285求調(diào)劑 +6	ytter 2026-03-22	6/300	2026-03-22 12:09 by 星空星月