collision.cpp

Lattice Boltzmann用于模拟方腔流的程序。 生成的plt文件需用tecplot软件打开。

#include "global.h"
void collision()
{
	register int x,y;
	
	register double rho,	/* density */
					vx,		/* x-velocity */						
					vy,		/* y-velocity */
					square,tau_inv,					
					f_eq0,f_eq1,f_eq2,f_eq3,f_eq4,f_eq5,f_eq6,f_eq7,f_eq8,/* 平衡态分布函数*/
					product;
					
	
	tau_inv=1.0/tau;
	
	for(y=0;y<max_y;y++)
	{
		for(x=0;x<max_x;x++)
		{
			if(wall[y][x]==FLUID)  /*流场内结点*/
			{
				rho=f[0][y][x]+f[1][y][x]+f[2][y][x]+f[3][y][x]+f[4][y][x]+f[5][y][x]+f[6][y][x]+f[7][y][x]+f[8][y][x]; //宏观量的计算
				vx=(f[1][y][x]-f[3][y][x]+f[5][y][x]+f[8][y][x]-f[6][y][x]-f[7][y][x])/rho;
				vy=(f[2][y][x]-f[4][y][x]+f[5][y][x]+f[6][y][x]-f[7][y][x]-f[8][y][x])/rho;
				
				/* 计算平衡态分布函数*/

				square =1.5*(vx * vx +vy *vy );
				f_eq0 =4./9.*rho*(1. - square);
	
				rho*=0.1111111111111111111111;
				f_eq1=rho*(1. + 3.0*vx + 4.5 *vx*vx - square);
				f_eq3=f_eq1-6.0*vx*rho;
				f_eq2=rho*(1. + 3.0*vy + 4.5 *vy*vy - square);
				f_eq4=f_eq2-6.0*vy*rho;
	
				rho*=0.25;

				product=vx+vy;
				f_eq5=rho* (1. + 3.0*product + 4.5 *product*product - square);
				f_eq7=f_eq5-6.0*product*rho;
				product=-vx+vy;
				f_eq6=rho* (1. + 3.0*product + 4.5 *product*product - square);
				f_eq8=f_eq6-6.0*product*rho;

				/*根据平衡态分布函数计算各节点的单粒子分布函数*/
			
				f[0][y][x]+=(f_eq0-f[0][y][x])*tau_inv;
				f[1][y][x]+=(f_eq1-f[1][y][x])*tau_inv;
				f[2][y][x]+=(f_eq2-f[2][y][x])*tau_inv;
				f[3][y][x]+=(f_eq3-f[3][y][x])*tau_inv;
				f[4][y][x]+=(f_eq4-f[4][y][x])*tau_inv;
				f[5][y][x]+=(f_eq5-f[5][y][x])*tau_inv;
				f[6][y][x]+=(f_eq6-f[6][y][x])*tau_inv;
				f[7][y][x]+=(f_eq7-f[7][y][x])*tau_inv;
				f[8][y][x]+=(f_eq8-f[8][y][x])*tau_inv;
			}
			else if(wall[y][x]==SET_U)  /*速度边界节点*/
			{
				rho=f[0][y][x]+f[1][y][x]+f[2][y][x]+f[3][y][x]+f[4][y][x]+f[5][y][x]+f[6][y][x]+f[7][y][x]+f[8][y][x];
				vx=u_0;
				vy=0.0;
				square =1.5*(vx * vx +vy *vy );
				f[0][y][x] =4./9.*rho*(1. - square);
	
				rho*=0.1111111111111111111111;
				f[1][y][x]=rho*(1. + 3.0*vx + 4.5 *vx*vx - square);
				f[3][y][x]=f[1][y][x]-6.0*vx*rho;
				f[2][y][x]=rho*(1. + 3.0*vy + 4.5 *vy*vy - square);
				f[4][y][x]=f[2][y][x]-6.0*vy*rho;
	
				rho*=0.25;

				product=vx+vy;
				f[5][y][x]=rho* (1. + 3.0*product + 4.5 *product*product - square);
				f[7][y][x]=f[5][y][x]-6.0*product*rho;
				product=-vx+vy;
				f[6][y][x]=rho* (1. + 3.0*product + 4.5 *product*product - square);
				f[8][y][x]=f[6][y][x]-6.0*product*rho;
			}
			else /* 固壁边界上的节点*/
			{
				rho=f[0][y][x]+f[1][y][x]+f[2][y][x]+f[3][y][x]+f[4][y][x]+f[5][y][x]+f[6][y][x]+f[7][y][x]+f[8][y][x];
				f[0][y][x]=rho*4.0/9.0;
			    f[1][y][x]=f[2][y][x]=f[3][y][x]=f[4][y][x]=rho/9.0;
			    f[5][y][x]=f[6][y][x]=f[7][y][x]=f[8][y][x]=rho/36.0;
			}
		}
	}
}