d2q9.cpp

格子Boltzmann方法 格子Boltzmann方法是为了保留格子气自动机方法的优点

// D2Q9.cpp
#include "StdAfx.h"		// VC++6.0 要求添加此行

#include "D2Q9.h"
#include "MyArrays.h"
#include <stdio.h>
#include <math.h>

int D=2, Q=9, Nx, Ny, n;
char e[9][2], e1[9], **boundary, bounceback, filename[80]="管道流.lbm";
double ***f, ***f1, **rho, **u, **v, *f0;
double mu, Re, omega, dx, dt, t, c, c_s;
double Lx, Ly, L0, U0, V0, Rho0, rho_max;

int MyUpper(char *str)		// 临时使用的处理字符串函数
{	// 删除str中的空格并将小写英文字母转换成大写字母，并将等号转换成'\0'，返回其后面的下标值
	int i=0, j=0, k=0;
	for(i=j=0; str[i]; i++)
	{
		if(str[i]!=' ')
		{
			if(str[i]>='a' && str[i]<='z')
				str[j] = str[i] - 'a'+'A';
			else if(str[i]=='=')
			{
				str[j] = '\0';
				k = j + 1;
			}
			else
				str[j] = str[i];
			j++;
		}
	}
	str[j] = '\0';
	return k;
}

void Init(char *filename)
{
	int i, j;
	char str[80];
	FILE *fp;

	e[0][0] = 0;  e[0][1] = 0;
	e[1][0] = 1;  e[1][1] = 0;
	e[2][0] = 0;  e[2][1] = 1;
	e[3][0] =-1;  e[3][1] = 0;
	e[4][0] = 0;  e[4][1] =-1;
	e[5][0] = 1;  e[5][1] = 1;
	e[6][0] =-1;  e[6][1] = 1;
	e[7][0] =-1;  e[7][1] =-1;
	e[8][0] = 1;  e[8][1] =-1;
	e1[0] = 0;
	e1[1] = 3;	e1[2] = 4;	e1[3] = 1;	e1[4] = 2;
	e1[5] = 7;	e1[6] = 8;	e1[7] = 5;	e1[8] = 6;

	if((fp=fopen(filename, "rt"))==NULL) exit(1);

	mu = 1.0;
	L0 = -1.0;
	bounceback = 0;
	while(1)
	{
		fgets(str, 80, fp);
		i = MyUpper(str);
		if(strcmp(str, "D")==0)
			sscanf(str+i, "%d", &D);
		else if(strcmp(str, "Q")==0)
			sscanf(str+i, "%d", &Q);
		else if(strcmp(str, "DOMAINSIZE")==0)
			sscanf(str+i, "%lf,%lf", &Lx, &Ly);
		else if(strcmp(str, "L0")==0)
			sscanf(str+i, "%lf", &L0);
		else if(strcmp(str, "BOUNCEBACK")==0)
			sscanf(str+i, "%d", &bounceback);
		else if(strcmp(str, "MESHSIZE")==0)
			sscanf(str+i, "%d,%d", &Nx, &Ny);
		else if(strcmp(str, "MU")==0)
			sscanf(str+i, "%lf", &mu);
		else if(strcmp(str, "RE")==0)
			sscanf(str+i, "%lf", &Re);
		else if(strcmp(str, "OMEGA")==0)
			sscanf(str+i, "%lf", &omega);
		else if(strcmp(str, "BOUNDARYDATA")==0)
			break;
	}
	if(L0<0) L0 = Ly;

	f   = (double***)Calloc3DArray(Nx, Ny, Q, sizeof(double));
	f1  = (double***)Calloc3DArray(Nx, Ny, Q, sizeof(double));
	rho = (double **)Calloc2DArray(Nx, Ny, sizeof(double));
	u   = (double **)Calloc2DArray(Nx, Ny, sizeof(double));
	v   = (double **)Calloc2DArray(Nx, Ny, sizeof(double));
	f0  = (double  *)Calloc1DArray(Q, sizeof(double));
	boundary=(char**)Calloc2DArray(Nx, Ny, sizeof(char));

	for(j=0; j<Ny; j++)
	{
		for(i=0; i<Nx; i++)		// 注意循环的顺序
		{
			fscanf(fp, "%c", str);
			boundary[i][j] = str[0] - '0';	
		}
		fgets(str, 80, fp);
	}
	fclose(fp);

	dx = Lx/Nx;
	dt = (2/omega-1)*dx*dx/(6*mu);
	c = dx/dt;
	c_s = c/sqrt(3);
	U0 = mu*Re/L0;
	V0 = 0;
	Rho0 = 1;

	t = 0;
	n = 0;
	init_macro();
	init_micro();
}

void Finish()
{
	if(boundary) Free2DArray((void**)boundary);
	if(f0)  Free1DArray((void*)f0);
	if(u)   Free2DArray((void**)v);
	if(v)   Free2DArray((void**)u);
	if(rho) Free2DArray((void**)rho);
	if(f1)  Free3DArray((void***)f1);
	if(f)   Free3DArray((void***)f);
}

void SetParameters()	// 仅改变参数 re, mu, omega 时，从头开始计算
{
	dt = (2/omega-1)*dx*dx/6.0/mu;
	c = dx/dt;
	c_s = c/sqrt(3);
	U0 = mu*Re/L0;
	V0 = 0;
	Rho0 = 1;

	t = 0;
	n = 0;
	init_macro();
	init_micro();
}

void init_macro()
{
	int i, j;
	for(i=0; i<Nx; i++)
	{
		for(j=0; j<Ny; j++)
		{
			if(boundary[i][j]==1)	// 流场流体入口处
			{
				u[i][j] = U0;
				v[i][j] = V0;
			}
			else					// 其它
				u[i][j] = v[i][j] = 0;
			rho[i][j] = Rho0;
		}
	}
}

void init_micro()
{
	int i, j, k;
	for(i=0; i<Nx; i++)
	{
		for(j=0; j<Ny; j++)
		{
			feq(u[i][j], v[i][j], rho[i][j], f0);
			for(k=0; k<Q; k++)
				f[i][j][k] = f0[k];
		}
	}
}

void feq(double u, double v, double rho, double *f0)
{
/*
	double u2, eu, eu2;
	int k;
	u /= c;
	v /= c;
	u2 = 1-1.5*(u*u + v*v);

	f0[0] = 4.0/9*rho*u2;
	for(k=1; k<Q; k++)
	{
		eu = e[k][0]*u + e[k][1]*v;
		eu2 = eu*eu;
		f0[k] = rho/9*(3*eu + 4.5*eu2 + u2);
	}
	for(k=5; k<Q; k++)
		f0[k] *= 0.25;
*/
	double temp_u2, temp_rho, uav, vmu;

	u /= c;
	v /= c;
	temp_u2 = 1-1.5*(u*u + v*v);
	temp_rho = 4.0*rho/9.0;
	uav = u+v;
	vmu = v-u;

	f0[0] = temp_rho * temp_u2;

	temp_rho *= 0.25;
	f0[1] = temp_rho * (temp_u2 + 3*u + 4.5*u*u);
	f0[2] = temp_rho * (temp_u2 + 3*v + 4.5*v*v);
	f0[3] = f0[1] - temp_rho * 6*u;
	f0[4] = f0[2] - temp_rho * 6*v;

	temp_rho *= 0.25;
	f0[5] = temp_rho * (temp_u2 + 3*uav + 4.5*uav*uav);
	f0[6] = temp_rho * (temp_u2 + 3*vmu + 4.5*vmu*vmu);
	f0[7] = f0[5] - temp_rho * 6*uav;
	f0[8] = f0[6] - temp_rho * 6*vmu;
}

void collision()
{
    int i, j, k;

	rho_max = 0;
    for(i=0; i<Nx; i++)
    {
        for(j=0; j<Ny; j++)
        {
			if(boundary[i][j]==1)
			{
				u[i][j] = U0;
				v[i][j] = V0;
				rho[i][j] = Rho0;
				continue;
			}
			else if(boundary[i][j]==3)
			{
				u[i][j] = 0;
				v[i][j] = 0;
				continue;
			}

            rho[i][j] = f[i][j][0];
			u[i][j] = v[i][j] = 0;
            for(k=1; k<Q; k++)
            {
                u[i][j] += e[k][0] * f[i][j][k];
                v[i][j] += e[k][1] * f[i][j][k];
                rho[i][j] += f[i][j][k];
            }
			u[i][j] *= c;
			v[i][j] *= c;
			if(rho[i][j]>rho_max) rho_max = rho[i][j];
			if(rho[i][j]>1e-6)
			{
				u[i][j] /= rho[i][j];
				v[i][j] /= rho[i][j];
			}
			else
			{
				u[i][j] = 0;
				v[i][j] = 0;
			}
		}
    }						// 计算宏观量

    for(i=0; i<Nx; i++)
    {
        for(j=0; j<Ny; j++)
        {
			feq(u[i][j], v[i][j], rho[i][j], f0);
            for(k=0; k<Q; k++)
                f1[i][j][k] = (1-omega)*f[i][j][k] + omega*f0[k];
        }
    }
}

void flow()
{
	int i, j, k, i1, j1;
	double f00[9];

	if(bounceback)	// 流动步，反弹边界处理
	{
		for(i=0; i<Nx; i++)
		{
			for(j=0; j<Ny; j++)
			{
				switch(boundary[i][j])
				{
				case 1:
					feq(u[i][j], v[i][j], rho[i][j], f1[i][j]); break;
				case 2:
					for(k=0; k<Q; k++)
						f1[i][j][k] = f1[i-1][j][k];			break;
				case 3:
					for(k=0; k<Q; k++)
					{
						i1 = i - e[k][0];
						j1 = j - e[k][1];
						if(i1>=0 && i1<Nx && j1>=0 && j1<Ny && boundary[i1][j1]==0)
							f1[i][j][k] = f1[i1][j1][e1[k]];
					}											break;
				}
			}
		}

		for(i=0; i<Nx; i++)
		{
			for(j=0; j<Ny; j++)
			{
				if(boundary[i][j]==0)
				{
					for(k=0; k<Q; k++)
						f[i][j][k] = f1[i-e[k][0]][j-e[k][1]][k];
				}
			}
		}
	}
	else		// 流动步，外推边界处理
	{
		for(i=0; i<Nx; i++)
	    {
		    for(j=0; j<Ny; j++)
			{
				switch(boundary[i][j])
				{
				case 0:
					for(k=0; k<Q; k++)
					{
						i1 = i - e[k][0];
						j1 = j - e[k][1];
						if(i1<0 || i1>=Nx || j1<0 || j1>=Ny)
						{
							i1 = i + e[k][0];
							j1 = j + e[k][1];
							if(i1<0 || i1>=Nx || j1<0 || j1>=Ny)	// 角点
							{
								f[i][j][k] = f1[i][j][k];
							}
							else	// 边点（外推）
							{
								feq(u[i][j], v[i][j], rho[i][j], f0);
								feq(u[i1][j1], v[i1][j1], rho[i1][j1], f00);
								f[i][j][k] = f1[i1][j1][k] + f0[k] - f00[k];
							}
						}
						else
							f[i][j][k] = f1[i1][j1][k];
					}
					break;
				case 1:
					feq(u[i][j], v[i][j], rho[i][j], f0);
					for(k=0; k<Q; k++)
						f[i][j][k] = f0[k];
					break;
				case 2:
					for(k=0; k<Q; k++)
						f[i][j][k] = f1[i-1][j][k];
					break;                        // 出口处理:无穷远边界
				default:
					feq(u[i][j], v[i][j], rho[i][j], f0);
					for(k=0; k<Q; k++)
						f[i][j][k] = f0[k];
					break;
				}
			}
		}
	}
}

void evolution()
{
	collision();
	flow();
	t += dt;
	n++;
}

void SaveFile()	// 将速度(u,v)、密度 rho 写入相应的文件
{
	FILE *fpu, *fpv, *fprho;
	int nx=Nx, ny=Ny;
	int i, j;

	fpu = fopen("u.dat", "wb");
	fpv = fopen("v.dat", "wb");
	fprho = fopen("rho.dat", "wb");

	fwrite(&nx, 1, sizeof(int), fpu);	// 按二进制格式存放（前2个整数皆为网格规模）
	fwrite(&ny, 1, sizeof(int), fpu);
	fwrite(u[0], nx*ny, sizeof(u[0][0]), fpu);

	fwrite(&nx, 1, sizeof(int), fpv);
	fwrite(&ny, 1, sizeof(int), fpv);
	fwrite(v[0], nx*ny, sizeof(v[0][0]), fpv);

	fwrite(&nx, 1, sizeof(int), fprho);
	fwrite(&ny, 1, sizeof(int), fprho);
	fwrite(rho[0], nx*ny, sizeof(rho[0][0]), fprho);

	fclose(fpu); fclose(fpv); fclose(fprho);

	fpu = fopen("u.txt", "wb");
	fpv = fopen("v.txt", "wb");
	fprho = fopen("rho.txt", "wb");

	for(j=0; j<Ny; j++)		// 写成文本文件格式
	{
		for(i=0; i<Nx; i++)
		{
			fprintf(fpu, "%f ", u[i][j]);
			fprintf(fpv, "%f ", v[i][j]);
			fprintf(fprho, "%f ", rho[i][j]);
		}
		fprintf(fpu, "\r\n");
		fprintf(fpv, "\r\n");
		fprintf(fprho, "\r\n");
	}
	fclose(fpu); fclose(fpv); fclose(fprho);	
}