// D2Q9.c[pp]#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, "\n");		fprintf(fpv, "\n");		fprintf(fprho, "\n");	}	fclose(fpu); fclose(fpv); fclose(fprho);	}