lbm模拟液滴从壁面滑落

#include
#include
#include
#include


#define R                                  18.0//initial（产生半圆液滴）
#define X0                                 1
#define Y0                                 （LY-30）

#define tau0                               1.0//流体性质参数
#define tau1                               1.0
#define m0                                 1.0
#define m1                                 1.0
#define G                                  0.2//作用力参数
#define GW0                                -0.01
#define GW1                                0.01
#define GR                                 -0.0004

//define
#define LX                                32
#define LY                                302
#define Q                                 9

#define fluid                              0//initial_geo
#define solid                              1



#define n0_in                              2.0
#define n0_out                             0.01
#define n1_in                              0.01
#define n1_out                             2.0



#define PI                                 3.1415926

#define ERR                                1.0e-6//收敛判断
#define IT                                 20000
#define max（xy）	                      （（（x） <= （y）） ? （y） : （x））

int flag[LX][LY];
double s[LX][LY] e[Q][2]
        f0[LX][LY][Q]f1[LX][LY][Q]
        n0[LX][LY]n1[LX][LY]n[LX][LY]den0[LX][LY]den1[LX][LY]den[LX][LY]
		veq0x[LX][LY]veq0y[LX][LY]veq1x[LX][LY]veq1y[LX][LY]vx[LX][LY]vy[LX][LY]//平衡态速度
		g0[LX][LY][Q]g1[LX][LY][Q]
		F0x[LX][LY]F0y[LX][LY]F1x[LX][LY]F1y[LX][LY]
	    pressure[LX][LY]p[LX][LY]
		C[LX][LY];
double rfa[2]gravity0gravity1;

double eq（int indexint kdouble ndouble uxdouble uy）;
void ini_geo（void）;
void initial（void）;
void collision（void）;
void stream（void）;
void para（void）;
void evolve（void）;
void contact_angle（int Tdouble errint N）;
//double err（double px[LX][LY]double py[LX][LY]）;
double err（double p[LX][LY]）;
void write（int T）;


double eq（int indexint kdouble ndouble uxdouble uy）
{
	double rf;
	rf=rfa[index];
	double euuueq;
	eu=ux*e[k][0]+uy*e[k][1];
	uu=ux*ux+uy*uy;
	if（k==0）
		eq=rf*n-2.0/3*n*uu;
	else if（k>=1&&k<=4）
		eq=（1-rf）*n/5+1.0/3*n*eu+1.0/2*n*eu*eu-1.0/6*n*uu;
	else
		eq=（1-rf）*n/20+1.0/12*n*eu+1.0/8*n*eu*eu-1.0/24*n*uu;

	return eq;
}
//ini_geo
void ini_geo（）
{
	int ij;
	for（j=0;j		for（i=0;i		{
			if（i==0||i==LX-1）
			{
				flag[i][j]=solid;
				s[i][j]=1.0;
			}
			else
			{
				flag[i][j]=fluid;
				s[i][j]=0.0;
			}
		}
}
//initial
void initial（）
{
	int ijk;
	for（j=0;j		for（i=0;i		{
			if（i==0||i==LX-1）
				n0[i][j]=n1[i][j]=0;
			else
			{
				if（（double（（i-X0）*（i-X0）+（j-Y0）*（j-Y0）））<=R*R）
				{
					n0[i][j]=n0_in;
					n1[i][j]=n1_in;
				}
				else
				{
					n0[i][j]=n0