//////////////////////////////////////////
// LBM MCMP Shan-Chen code //
//////////////////////////////////////////
// Dong-Ke SUN //
// dongke.sun@gmail.com //
// www.calcflows.com //
// 2012.07.07 //
//////////////////////////////////////////
#ifndef CLBMCMPSC_H
#define CLBMCMPSC_H
#include"LBMStudio.h"
CLBMcmpSC::CLBMcmpSC()
{
m_ux=Define2DArray(Nx,Ny);
m_uy=Define2DArray(Nx,Ny);
F_x=Define2DArray(Nx,Ny);
F_y=Define2DArray(Nx,Ny);
G_x=Define2DArray(Nx,Ny);
G_y=Define2DArray(Nx,Ny);
Fx_c=Define2DArray(Nx,Ny);
Fy_c=Define2DArray(Nx,Ny);
psi_c=Define2DArray(Nx,Ny);
phi=Define2DArray(Nx,Ny);
bubble= Define2DintArray(Nx,Ny);
fclose(denfp);
denfp = fopen("density.txt","wt");
fprintf(denfp,"frame,\taveDensity,\tavePhi,\taveG,\tr_ave,\tr_max,\tr_min,\tn_bubble*\n");
InitialGeography();
m_rho_f_0 = m_rho_f;
m_rho_g_0 = m_rho_g;
density_ratio = m_rho_f/m_rho_g;
m_rho_f /= m_rho_f_0;
m_rho_g /= m_rho_g_0;
//{
// printf("density_ratio=%g\n",density_ratio);
// system("pause");
//}
char pdfType='R';
switch(pdfType)
{
case 'R':
SetRandPDF();
break;
case 'U':
SetUniformPDF();
// GenerateOneSeed();
break;
default:
SetRandPDF();
;
}
// GenerateOneSeed(); // for only one bubble
//WriteResult(1000);
// GenerateSeeds(100); // for multi bubbles
}
CLBMcmpSC::~CLBMcmpSC()
{
Delete2DArray(m_ux,Nx);
Delete2DArray(m_uy,Nx);
Delete2DArray(F_x,Nx);
Delete2DArray(F_y,Nx);
Delete2DArray(G_x,Nx);
Delete2DArray(G_y,Nx);
Delete2DArray(Fx_c,Nx);
Delete2DArray(Fy_c,Nx);
Delete2DArray(psi_c,Nx);
Delete2DArray(phi,Nx);
Delete2DArray(bubble,Nx);
}
void CLBMcmpSC::GenerateOneSeed( )
{
int x, y;
//for single bubble
int R=Nx/4;//5;//4;//
int R_squ=R*R;//1;//
double xCenter=Nx/2;//
double yCenter=Ny/2;//
double num1;
double num2;
for(x=0;x<Nx;x++)
{
for(y=0; y<Ny; y++)
{
//if one bubble growth is to be studied, switch the following two lines
num1 = 0.05;// (num1+1.0)/2.0;//= (rand())/(1.0*RAND_MAX);//suiji fenbu
num2 = 0.95;//num2/2.0;//
if((x-xCenter)*(x-xCenter)+(y-yCenter)*(y-yCenter)<R_squ)
{
CalEquStandard(f_out[x][y],0,0,num1*m_rho_f);
CalEquStandard(g_out[x][y],0,0,num2*m_rho_g);
//CalEquSC(f_out[x][y],0,0,0.05*m_rho_f);
//CalEquSC(g_out[x][y],0,0,0.95*m_rho_f);
phi[x][y]=0.025;
}
else
{
//CalEquSC(f_out[x][y],0,0,m_rho_g);
//CalEquSC(g_out[x][y],0,0,m_rho_f);
}
/* */
/*
if((x-xCenter)*(x-xCenter)+(y-yCenter)*(y-yCenter)<R_squ)
{
if(!is_obst[x][y])
{
CalEquStandard(f_out[x][y],0,0,0.05*m_rho_f);
CalEquStandard(g_out[x][y],0,0,0.95*m_rho_f);
phi[x][y]=0.025;
}
}
*/
/*
if((x-2*Nx/4)*(x-2*Nx/4)+(y-Ny/2)*(y-Ny/2)<R_squ)
{
CalEquStandard(f_out[x][y],0,0,0.05*m_rho_f);
CalEquStandard(g_out[x][y],0,0,0.95*m_rho_f);
phi[x][y]=0.025;
}
if((x-3*Nx/4)*(x-3*Nx/4)+(y-Ny/2)*(y-Ny/2)<R_squ)
{
CalEquStandard(f_out[x][y],0,0,0.05*m_rho_f);
CalEquStandard(g_out[x][y],0,0,0.95*m_rho_f);
}*/
}
}
}
void CLBMcmpSC::GenerateSeeds(int total_bubble)
{
int x, y;
int bx, by; // bubble node location
int bx0,by0;// bubble center location
int R=2;
int R_squ=R*R;
//检查total_bubble设置是否合理
{
int sum=Nx*Ny;
if(4*R*R*total_bubble >= sum/3)
{
printf("Attention!!!\n");
printf("Total number of bubbles is not reasonable!\n");
exit(0);
}
}
srand( (unsigned)time(NULL) ); //srand()函数产生一个以当前时间开始的随机种子
for(x=0; x<Nx; x++)
{
for(y=0; y<Ny; y++)
{
CalEquStandard(f_out[x][y],0,0,0.95*m_rho_f);
CalEquStandard(g_out[x][y],0,0,0.05*m_rho_f);
if(!is_obst[x][y])
{
CalEquStandard(c_out[x][y],0,0,m_C0);
}
else
{
CalEquStandard(c_out[x][y],0,0,0);
}
phi[x][y]=1.0;
bubble[x][y]=0;
}
}
int counter = 0;//计数
srand( (unsigned)time(NULL) ); //srand()函数产生一个以当前时间开始的随机种子
while(counter<total_bubble)
{
int is_OK=FALSE; //检查随机产生的bubble圆心是否合适
while(!is_OK)
{
//随机产生一点坐标
bx0=Nx*(rand())/(RAND_MAX);
by0=Ny*(rand())/(RAND_MAX);
is_OK=TRUE;
//检查以此点为圆心所囊括范围内所有的点是否满足
//1. 不是obst
//2. 不是bubble
for(bx=-R; bx<=R; bx++)
{
x = bx0+bx;
for(by=-R; by<=R; by++)
{
y=by0+by;
if((x-bx0)*(x-bx0)+(y-by0)*(y-by0)<=R_squ+2)
{
if(x<0) x = x+Nx;
if(x>=Nx)x = x%Nx;
if(y<0) y = y+Ny;
if(y>=Ny)y = y%Ny;
if(bubble[x][y]==1||is_obst[x][y])
{
is_OK=FALSE;
break;
}
}
}
if(!is_OK)
{
break;
}
}
}
if(is_OK)
{
counter++;
for(bx=-R; bx<=R; bx++)
{
x=bx0+bx;
for(by=-R; by<=R; by++)
{
y=by0+by;
if((x-bx0)*(x-bx0)+(y-by0)*(y-by0)<R_squ)
{
if(x<0) x = x+Nx;
if(x>=Nx)x = x%Nx;
if(y<0) y = y+Ny;
if(y>=Ny)y = y%Ny;
bubble[x][y]=1;
double f_loc=0.05*m_rho_f;
double g_loc=0.95*m_rho_f;
CalEquStandard(f_out[x][y],0,0,f_loc);
CalEquStandard(g_out[x][y],0,0,g_loc);
phi[x][y] = (f_loc-g_loc)/(g_loc+f_loc)/2.0 + 0.5;
// EquilibriumF(c_out[x][y],0,0,0.22);
}
}
}
}
}
}
void CLBMcmpSC::InitialGeography()
{
int x, y;
//for single bubble
int R=Ny/5;//4;//
int R_squ=R*R;//1;//
double xCenter=Nx/2;
double yCenter=Ny/2;//Ny-R-5;
for(x=0; x<Nx; x++)
{
for(y=0; y<Ny; y++)
{
m_solidFrc[x][y]=0;
}
}
// //环形流道
// double x_o=Nx/2;
// double y_o=Ny/2;
// double r=Nx/4;
////#pragma omp parallel for private(i,x,y)
// for(x=0;x<Nx;x++)
// {
// for(y=0;y<Ny;y++)
// {
// //double d = dist(x,y_o, x,y);
// //if(d >= Ny/2.0-3)
// double d = dist(x_o,y_o, x,y);
// if(d >= Nx/2.0 || d<=r)
// {
// m_solidFrc[x][y] = 1;
// }
// else
// {
// m_solidFrc[x][y] = 0;
// }
//
// }
// }
// //螺旋流道
// int x_o=Nx/2;
// int y_o=Ny/2;
// double r=10;
//#pragma omp parallel for private(i,x,y)
// for(x=0;x<Nx;x++)
// {
// for(y=0;y<Ny;y++)
// {
// double t(0);
// double b=2;
// while(t<100)
// {
// double x_a = x_o + b*t*cos(1.0*t);
// double y_a = y_o + b*t*sin(1.0*t);
// double d = dist(x_a,y_a, x,y);
//
// if(d <= 1.5)
// {
// m_solidFrc[x][y] = 1;
// break;
// }
// else
// {
// m_solidFrc[x][y] = 0;
// }
//
// t += 0.02;
// }
// }
// }
SetFlags();
}
double ** CLBMcmpSC::GetPhiField()
{
return phi;
}
double** CLBMcmpSC::GetCField()
{
int x,y;
for(x=0; x<Nx; x++)
{
for(y=0; y<Ny; y++)
{
m_C_new[x][y]=GetDensity(c_out[x][y],QC);
}
}
return m_C_new;
}
void CLBMcmpSC::Bounceback()
{
int x,y;
for(x = 0;x< Nx ; x++)
{
for(y = 0;y< Ny;y++)
{
//if (is_obst[x][y])
if(is_obst[x][y]||is_inter[x][y])//如果遇到固相或者障碍物,分布函数沿原路返回
{
f_out[x][y][1] = f_hlp[x][y][3];//c...........east
f_out[x][y][2] = f_hlp[x][y][4];//c...........north
f_out[x][y][3] = f_hlp[x][y][1];//c...........west
f_out[x][y][4] = f_hlp[x][y][2];//c...........south
f_out[x][y][5] = f_hlp[x][y][7];//c...........north-east
f_out[x][y][6] = f_hlp[x][y][8];//c...........north-west
f_out[x][y][7] = f_hlp[x][y][5];//c...........south-west
f_out[x][y][8] = f_hlp[x][y][6];//c..
ShanChen_D2Q9_MCMP.zip_LBM_lbm shanchen_mcmp_圆柱绕流_绕流
版权申诉
5星 · 超过95%的资源 10 浏览量
2022-07-15
16:46:04
上传
评论
收藏 1.71MB ZIP 举报
寒泊
- 粉丝: 75
- 资源: 1万+
最新资源
- 基于Matlab人脸肤色定理的教师人数统计+源代码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
- 基于Matlab霍夫曼变换的表盘读数识别+源代码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
- 基于Matlab火灾烟雾检测源码带GUI界面+源代码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
- 基于Matlab的恶劣天气交通标志识别系统+源代码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
- 基于MATLAB的霍夫曼变换的表盘示数识别+源代码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
- 基于Matlab的车道线识别系统 +源代码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
- 基于MATLAB的教室人数统计系统带Gui界面+源代码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
- 基于MATLAB的教室人数统计系统带Gui界面+源代码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
- 基于MATLAB 的霍夫曼变换答题卡识别源码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
- 基于Matlab+bp神经网络的神经网络汉字识别系统+源代码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
评论3