about Spurious currents in SCMP. (with C++ code)

Lattice Boltzmann Method LBM: Theory
#1

Hi,

Recnely, I study the SCMP at large density ratio with lattice Boltzmann method. My study base on Peng Yuan and Laura Schaefer’s work: Equations of state in a lattice Boltzmann model. In this paper, authors research different equations of state with spurious currents and the density ratio. The spurious current is smaller than 0.1 in most cases, except spurious currents in large density ratio case(>1000) and vdW EOS case is larger than 0.1. In my simulation, spurious currents is unusually large(>=0.2) in small density ratio or large density ratio. For example, when I set temperature T=0.06, the density ratio is about 195 and the spurious currents maximun is 0.21076, when I set temperature T=0.07, the density ratio is about 34 and the spurious currents maximun is 0.133204. Compared to the results in paper, the deviation is large.

I recheck my code and model many times. I don’t know this phenomenon is caused by my fault or the model? Anyone use this model or encounter this phenomenon? And I attach my code on the topic.

Thanks for your time.

Kind regards,
Skying

#2

//SCMP with large density ration

#include
#include
#include
#include <stdio.h>
using namespace std;

//Parameters
const int num=5e4;	//total number of iterations
const int xnum=128;	
const int ynum=128;
const int x=64;
const int y=64;
const int r=12;
const double cs=pow(3.,-0.5);	//lattice sound speed
const double r_t=1.;	//relax time
const double g=-1.;		//interparticle force parameter
const double tim=1.;	//time step
const double tem=0.07;		//temperature
const double a=1.;			// C-S EOS parameters
const double b=4.;
const double R=1.;
const double denseL=0.33768;		//initial density of liquid phase
const double denseV=1.05e-2;	//initial density of vapor phase
const int ex[9]={0,1,0,-1,0,1,-1,-1,1};	
const int ey[9]={0,0,1,0,-1,1,1,-1,-1};

double dist(double dense,double ux,double uy,int ex,int ey,int mode); //calculate distribution funtions
void text(double ux[xnum][ynum],double uy[xnum][ynum],double dense[xnum][ynum],double press[xnum][ynum],double velocity[xnum][ynum],int t); //data output

int main(){

int i,j,in,jn,mode;					//coordinate
int t;	
double ux[xnum][ynum];				//equilibrium velocity in X direction
double uy[xnum][ynum];				//equilibrium velocity in Y direction
double velocity[xnum][ynum];
double dense[xnum][ynum];	
double press[xnum][ynum];	
double psi[xnum][ynum];				//effective mass
double fx[xnum][ynum];				//force in X direction
double fy[xnum][ynum];				//force in X direction
double dis[xnum][ynum][9];			//distribution funtions
double colli[xnum][ynum][9];		//equilibrium distribution funtions
double stream[xnum][ynum][9];		//distribution funtions streaming

//density initialization
for(i=0;i<xnum;i++){
	for(j=0;j<ynum;j++){
		//dense[i][j]=(rand()%1000)/1e6;
		//dense[i][j]+=0.1;

		if((pow(i-x,2.)+pow(j-y,2.))<=r*r)
			dense[i][j]=denseL;
		else
			dense[i][j]=denseV;

		ux[i][j]=0;
		uy[i][j]=0;
		for(mode=0;mode<9;mode++)
			dis[i][j][mode]=dist(dense[i][j],ux[i][j],uy[i][j],ex[mode],ey[mode],mode);
	}
}

//interation
for(t=0;t<=num;t++){

	//streaming
	for(i=0;i<xnum;i++){
		for(j=0;j<ynum;j++){
			for(mode=0;mode<9;mode++){
				if(i+ex[mode]<0)
					in=xnum-1;
				else if(i+ex[mode]>=xnum)
					in=0;
				else
					in=i+ex[mode];

				if(j+ey[mode]<0)
					jn=ynum-1;
				else if(j+ey[mode]>=ynum)
					jn=0;
				else
					jn=j+ey[mode];

				stream[in][jn][mode]=dis[i][j][mode];
			}
		}
	}

	//update dense, velocity, effective mass and pressure
	for(i=0;i<xnum;i++){
		for(j=0;j<ynum;j++){
			dense[i][j]=0.;
			ux[i][j]=0.;
			uy[i][j]=0.;
			for(mode=0;mode<9;mode++){
				dense[i][j]+=stream[i][j][mode];
				ux[i][j]+=ex[mode]*stream[i][j][mode];
				uy[i][j]+=ey[mode]*stream[i][j][mode];
			}
			ux[i][j]/=dense[i][j];
			uy[i][j]/=dense[i][j];

			press[i][j]=dense[i][j]*R*tem*(1+b*dense[i][j]/4.+pow(b*dense[i][j]/4.,2)-pow(b*dense[i][j]/4.,3))/pow(1-b*dense[i][j]/4.,3)-a*pow(dense[i][j],2);
			psi[i][j]=pow(2.*(press[i][j]-cs*cs*dense[i][j])/(6.*g),0.5);
			//psi[i][j]=pow(6.*(press[i][j]-cs*cs*dense[i][j])/g,0.5);
			//psi[i][j]=4.*exp(-200./dense[i][j]);
		}
	}

	//interparticle force calculation
	for(i=0;i<xnum;i++){
		for(j=0;j<ynum;j++){
			fx[i][j]=0;
			fy[i][j]=0;
			for(mode=1;mode<=4;mode++){
				if(i+ex[mode]<0)
					in=xnum-1;
				else if(i+ex[mode]>=xnum)
					in=0;
				else
					in=i+ex[mode];

				if(j+ey[mode]<0)
					jn=ynum-1;
				else if(j+ey[mode]>=ynum)
					jn=0;
				else
					jn=j+ey[mode];

				fx[i][j]+=ex[mode]*psi[in][jn]/3.;
				fy[i][j]+=ey[mode]*psi[in][jn]/3.;
			}

			for(mode=5;mode<=8;mode++){
				if(i+ex[mode]<0)
					in=xnum-1;
				else if(i+ex[mode]>=xnum)
					in=0;
				else
					in=i+ex[mode];

				if(j+ey[mode]<0)
					jn=ynum-1;
				else if(j+ey[mode]>=ynum)
					jn=0;
				else
					jn=j+ey[mode];

				fx[i][j]+=ex[mode]*psi[in][jn]/12.;
				fy[i][j]+=ey[mode]*psi[in][jn]/12.;
			}

			fx[i][j]*=-6.*g*psi[i][j];
			fy[i][j]*=-6.*g*psi[i][j];
		}
	}

	//collision
	for(i=0;i<xnum;i++){
		for(j=0;j<ynum;j++){
			for(mode=0;mode<9;mode++){
				colli[i][j][mode]=dist(dense[i][j],ux[i][j]+fx[i][j]*r_t/dense[i][j],uy[i][j]+fy[i][j]*r_t/dense[i][j],ex[mode],ey[mode],mode);
				dis[i][j][mode]=stream[i][j][mode]-(stream[i][j][mode]-colli[i][j][mode])/r_t;
			}
			velocity[i][j]=pow(ux[i][j]+0.5*fx[i][j]/dense[i][j],2.)+pow(uy[i][j]+0.5*fy[i][j]/dense[i][j],2.);
			velocity[i][j]=pow(velocity[i][j],0.5);
		}
	}

	if(fmod(t,200.)==0)
		cout<<t<<"\t"<<dense[xnum/2][ynum/2]<<"\n";
	if(fmod(t,10000.)==0.)
		text(ux,uy,dense,press,velocity,t);
}

//text(ux,uy,dense,press,velocity,t);

return 0;

}

double dist(double dense,double ux,double uy,int ex,int ey,int mode){
	double dis;

	if(mode==0)
		dis=4.*dense*(1-1.5*(ux*ux+uy*uy))/9.;
	else if(mode>4)
		dis=dense*(1.+3.*(ex*ux+ey*uy)+4.5*pow(ex*ux+ey*uy,2)-1.5*(ux*ux+uy*uy))/36.;
	else
		dis=dense*(1.+3.*(ex*ux+ey*uy)+4.5*pow(ex*ux+ey*uy,2)-1.5*(ux*ux+uy*uy))/9.;

	return dis;
}