Dear all ,
I am new to this forum, and want to get help to correct the code I have written,
This is very simple lid driven cavity code using D2Q9 model.
The Zou and He boundary condition is used at moving lid. The compilation of code is fine but after running the
code I am getting nan values, which is due to reduction in density at moving lid,
but I have coded according to Zou and He paper but still getting error.
Can anybody please help me to site the error in it.
Thanks in advance
Bug found!
code corrected.
Dear Atulbhagat,
You are probably dividing by 0 somewhere in the code.
But if you can show the code we might be able to fix it.
Regards
Githin
i am working on a project fdlbm(finite difference latiice boltezmam )
and i need a code for lid cavity
and i need it very much
best regards
mahmoud
Hi mahmoud,
This is the code I have written rather say converted from fortran to c. This code you can get in A.A.Mohammad Book.
Hope this helps
[code=“cpp”]
#include
#include
#include
#include
#include <time.h>
const int nx = 100, ny = 100, q = 9;
double csq, w0 = 4.0/9.0, w1 = 1.0/9.0, w2 = 1.0/36.0;
double ***f, ***feq, **rho, **u, **v, **strf;
double dx, dt, alpha, omega, rhon, u0 = 0.1, v0 = 0.0;
double R, U, V, u2, v2, uv, s2, uv2, uvm, uvm2;
FILE *data2d;
// Memory allocation
void memoryallocation()
{
feq = new double **[q];
f = new double **[q];
for (int i = 0; i < q; i++)
{
feq[i] = new double *[nx];
f[i] = new double *[nx];
for (int j = 0; j < nx; j++)
{
feq[i][j] = new double [ny];
memset (feq[i][j], 0.0, (ny)*sizeof(double));
f[i][j] = new double [ny];
memset (f[i][j], 0.0, (ny)*sizeof(double));
}
}
rho = new double *[nx];
u = new double *[nx];
v = new double *[nx];
strf = new double *[nx];
for (int i = 0; i < nx; i++)
{
rho[i] = new double [ny];
memset (rho[i], 0.0, (ny)*sizeof(double));
u[i] = new double [ny];
memset (u[i], 0.0, (ny)*sizeof(double));
v[i] = new double [ny];
memset (v[i], 0.0, (ny)*sizeof(double));
strf[i] = new double [ny];
memset (strf[i], 0.0, (ny)*sizeof(double));
}
}
// Deletion of Memory
void deletememory()
{
for (int i = 0; i < q; i++)
{
for (int j = 0; j < nx; j++)
{
delete[] feq[i][j];
delete[] f[i][j];
}
delete[] feq[i];
delete[] f[i];
}
for (int i = 0; i < nx; i++)
{
delete[] rho[i];
delete[] u[i];
delete[] v[i];
delete[] strf[i];
}
}
// Equilibrium state
void equilibrium()
{
for(int i = 0; i < nx; i++)
{
for(int j = 0; j < ny; j++)
{
U = u[i][j];
V = v[i][j];
R = rho[i][j];
u2 = UU;
v2 = VV;
s2 = -1.5*(u2+v2);
uv = U+V;
uv2 = (U+V)(U+V);
uvm = U-V;
uvm2 = (U-V)(U-V);
feq[0][i][j] = w0R(1. + s2);
feq[1][i][j] = w1R(1. + 3.U + 4.5u2 + s2);
feq[2][i][j] = w1R(1. + 3.V + 4.5v2 + s2);
feq[3][i][j] = w1R(1. - 3.U + 4.5u2 + s2);
feq[4][i][j] = w1R(1. - 3.V + 4.5v2 + s2);
feq[5][i][j] = w2R(1. + 3.uv + 4.5uv2 + s2);
feq[6][i][j] = w2R(1. - 3.uvm + 4.5uvm2 + s2);
feq[7][i][j] = w2R(1. - 3.uv + 4.5uv2 + s2);
feq[8][i][j] = w2R(1. + 3.uvm + 4.5uvm2 + s2);
}
}
}
// Initial conditions
void initialstate()
{
for (int i = 0; i < nx; i++)
{
for (int j = 0; j < ny; j++)
{
rho[i][j] = 5.0;
strf[i][j] = 0.0;
}
}
equilibrium();
for(int k = 0; k < q; k++)
{
for(int i = 0; i < nx; i++)
{
for(int j = 0; j < ny; j++)
{
f[k][i][j] = feq[k][i][j];
}
}
}
}
// Macro-variable calculations
void macrovar()
{
for(int i = 0; i < nx; i++)
{
for(int j = 0; j < ny; j++)
{
rho[i][j] = 0.0;
for(int k = 0; k < q; k++)
{
rho[i][j] += f[k][i][j];
}
}
}
for(int i = 0; i < nx; i++)
{
for(int j = 0; j < ny; j++)
{
u[i][j] = (f[1][i][j] - f[3][i][j] + f[5][i][j] - f[6][i][j] - f[7][i][j] + f[8][i][j])/rho[i][j];
v[i][j] = (f[5][i][j] + f[2][i][j] + f[6][i][j] -f[7][i][j] -f[4][i][j] - f[8][i][j])/rho[i][j];
}
}
}
// collision
void collision()
{
for(int k = 0; k < q; k++)
{
for(int i = 0; i < nx; i++)
{
for(int j = 0; j < ny; j++)
{
f[k][i][j] = f[k][i][j](1.0-omega) + omegafeq[k][i][j];
}
}
}
}
//streaming
void streaming()
{
for(int i = (nx-1); i > 0; i–)
{
for(int j = (ny-1); j > 0; j–)
{
f[1][i][j] = f[1][i-1][j];
f[5][i][j] = f[5][i-1][j-1];
}
}
for(int i = 0; i < (nx-1); i++)
{
for(int j = 0; j < (ny-1); j++)
{
f[3][i][j] = f[3][i+1][j];
f[7][i][j] = f[7][i+1][j+1];
}
}
for(int j = (ny-1); j > 0; j–)
{
for(int i = 0; i < (nx-1); i++)
{
f[2][i][j] = f[2][i][j-1];
f[6][i][j] = f[6][i+1][j-1];
}
}
for(int i = (nx-1); i > 0; i–)
{
for(int j = 0; j < (ny-1); j++)
{
f[4][i][j] = f[4][i][j+1];
f[8][i][j] = f[8][i-1][j+1];
}
}
}
// boundary condition
void boundarycondition()
{
for(int j = 0; j < ny; j++)
{
// left layer i.e. west boundary
f[1][0][j] = f[3][0][j];
f[5][0][j] = f[7][0][j];
f[8][0][j] = f[6][0][j];
// right layer i.e. east boundary
f[3][nx-1][j] = f[1][nx-1][j];
f[6][nx-1][j] = f[8][nx-1][j];
f[7][nx-1][j] = f[5][nx-1][j];
}
// bottom layer i.e. south boundary
for(int i = 0; i < nx; i++)
{
f[2][i][0] = f[4][i][0];
f[5][i][0] = f[7][i][0];
f[6][i][0] = f[8][i][0];
}
// top layer i.e. north boundary (a moving plate/lid)
for(int i = 1; i < (nx-1); i++)
{
rhon = (1.0/(1.0+v0))(f[0][i][ny-1]+f[1][i][ny-1]+f[3][i][ny-1]+2.0(f[2][i][ny-1]+f[6][i][ny-1]+f[5][i][ny-1]));
f[4][i][ny-1] = f[2][i][ny-1] - ((2.0/3.0)rhonv0);
f[7][i][ny-1] = f[5][i][ny-1]-(rhonv0/6.0)-((1.0/2.0)rhonu0);
f[8][i][ny-1] = f[6][i][ny-1]-(rhonv0/6.0)+((1.0/2.0)rhonu0);
}
}
int main()
{
// Time calculations
clock_t t;
t = clock();
// Variable declaration
dt = 1.0;
dx = 1.0;
alpha = 0.01;
// Parameter calculation
csq = dxdx/(dtdt);
omega = 1.0/(3.alpha/(dtcsq)+0.5);
memoryallocation();
initialstate();
int nstep = 200000;
// main loop of LBM
for (int m = 0; m < nstep; m++)
{
equilibrium();
collision();
streaming();
boundarycondition();
macrovar();
if(m % 100 == 0)
{
printf(“nsteps = %d , U(nx/2,ny/2) = %e\n”, m, u[nx/2][ny/2]);
}
}
data2d = fopen(“lid_mid_plane_output.txt”, “w”);
for (int i = 0; i < nx; i++)
{
fprintf(data2d, “%i\t%f\t%f\t%f\n”, i, rho[i][ny/2], u[i][ny/2], v[i][ny/2]); //X,y,rho,u,v;
}
fclose(data2d);
for (int i = 1; i < nx; i++)
{
double rhoav = 0.5*(rho[i-1][0]+rho[i][0]);
if (i != 0) strf[i][0] = strf[i-1][0] - rhoav0.5(v[i-1][0]+v[i][0]);
for (int j = 1; j < ny; j++)
{
double rhom = 0.5*(rho[i][j]+rho[i][j-1]);
strf[i][j] = strf[i][j-1] + rhom*0.5*(u[i][j-1]+u[i][j]);
}
}
data2d = fopen(“stremfun_output.txt”, “w”);
for(int j = 0; j < ny; j++)
{
for (int i = 0; i < nx; i++)
{
fprintf(data2d, “%f\t”, strf[i][j]);
}
fprintf(data2d, “\n”);
}
fclose(data2d);
deletememory();
t = clock() - t;
double time_taken = ((double)t)/CLOCKS_PER_SEC;
printf (“Time of execution = %f\n”, time_taken);
return 0;
}
hi AtulBhagat
sorry for replying late,thanks for the code
but the code i need is finite difference lbm code
which have not streaming in that method
once again very very thanks for the code
if you have a code for linearly heating natural convection cavity
i run the heating natural convection cavity at A.A.Mohammad Book. and it work fine
and when i modify it to run linearly heating natural convection cavity it did not work fine
and i do not why
i changed the boundary and i think i changed write
6.6.1 Example: Natural Convection in a Differentially Heated Cavity page 98 and it run with me good
and here is my code
! this code for linealy heated cavity with all wall at static with linealy left heated wall and adiabatic at north
! and heated at bottom and low temp at right
[code="fortran"parameter (n=100,m=100)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8), cx(0:8),cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
real g(0:8,0:n,0:m), geq(0:8,0:n,0:m),th(0:n,0:m)
integer i
open(2,file=‘uvfield’)
open(3,file=‘uvely’)
open(4,file=‘tmx’)
cx(:)=(/0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0/)
cy(:)=(/0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0/)
w(:)=(/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./)
uo=0.0
sumvelo=0.0
rhoo=6.00
dx=1.0
dy=dx
dt=1.0
tw=1.0
th=0.0
ra=1.0e3
pr=0.7
visco=0.02
alpha=visco/pr
pr=visco/alpha
gbeta=raviscoalpha/(float(mmm))
Re=uom/alpha
print , “Re=”, Re
omega=1.0/(3.visco+0.5)
omegat=1.0/(3.alpha+0.5)
mstep=100000
do j=0,m
do i=0,n
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
end do
do i=0,n
u(i,m)=uo
v(i,m)=0.0
end do
! main loop
1 do kk=1,mstep
call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m,th,gbeta)
call streaming(f,n,m)
call bounceb(f,n,m)
call rhouv(f,rho,u,v,cx,cy,n,m)
! ——————————–
! collestion for scalar
call collt(u,v,g,geq,th,omegat,w,cx,cy,n,m)
! streaming for scalar
call streaming(g,n,m)
call gbound(g,tw,w,n,m)
do j=0,m
do i=0,n
sumt=0.0
do k=0,8
sumt=sumt+g(k,i,j)
end do
th(i,j)=sumt
end do
end do
!print , th(n/2,m/2),v(5,m/2),rho(0,m/2),u(n/2,m/2),v(n/2,m/2),rho(n,m/2)
END DO
! end of the main loop
call result(u,v,rho,th,uo,n,m,ra)
stop
end
! end of the main program
subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m,th,gbeta)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8), cx(0:8),cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
real th(0:n,0:m)
tref=0.50
DO i=0,n
DO j=0,m
t1=u(i,j)u(i,j)+v(i,j)v(i,j)
DO k=0,8
t2=u(i,j)cx(k)+v(i,j)cy(k)
force=3.w(k)gbeta(th(i,j)-tref)cy(k)rho(i,j)
if(i.eq.0.or.i.eq.n) force =0.0
if(j.eq.0.or.j.eq.m) force =0.0
feq(k,i,j)=rho(i,j)w(k)(1.0+3.0t2+4.50t2t2-1.50t1)
f(k,i,j)=omegafeq(k,i,j)+(1.-omega)f(k,i,j)+force
END DO
END DO
END DO
return
end
subroutine collt(u,v,g,geq,th,omegat,w,cx,cy,n,m)
real g(0:8,0:n,0:m),geq(0:8,0:n,0:m),th(0:n,0:m)
real w(0:8),cx(0:8),cy(0:8)
real u(0:n,0:m),v(0:n,0:m)
do i=0,n
do j=0,m
do k=0,8
geq(k,i,j)=th(i,j)w(k)(1.0+3.0(u(i,j)cx(k)+v(i,j)cy(k)))
g(k,i,j)=omegatgeq(k,i,j)+(1.0-omegat)g(k,i,j)
end do
end do
end do
return
end
subroutine streaming(f,n,m)
real f(0:8,0:n,0:m)
! streaming
DO j=0,m
DO i=n,1,-1 !RIGHT TO LEFT
f(1,i,j)=f(1,i-1,j)
END DO
DO i=0,n-1 !LEFT TO RIGHT
f(3,i,j)=f(3,i+1,j)
END DO
END DO
DO j=m,1,-1 !TOP TO BOTTOM
DO i=0,n
f(2,i,j)=f(2,i,j-1)
END DO
DO i=n,1,-1
f(5,i,j)=f(5,i-1,j-1)
END DO
DO i=0,n-1
f(6,i,j)=f(6,i+1,j-1)
END DO
END DO
DO j=0,m-1 !BOTTOM TO TOP
DO i=0,n
f(4,i,j)=f(4,i,j+1)
END DO
DO i=0,n-1
f(7,i,j)=f(7,i+1,j+1)
END DO
DO i=n,1,-1
f(8,i,j)=f(8,i-1,j+1)
END DO
END DO
return
end
subroutine bounceb(f,n,m)
real f(0:8,0:n,0:m)
do j=0,m
!west boundary
f(1,0,j)=f(3,0,j)
f(5,0,j)=f(7,0,j)
f(8,0,j)=f(6,0,j)
!east boundary
f(3,n,j)=f(1,n,j)
f(7,n,j)=f(5,n,j)
f(6,n,j)=f(8,n,j)
end do
do i=0,n
!south boundary
f(2,i,0)=f(4,i,0)
f(5,i,0)=f(7,i,0)
f(6,i,0)=f(8,i,0)
!north boundary
f(4,i,m)=f(2,i,m)
f(8,i,m)=f(6,i,m)
f(7,i,m)=f(5,i,m)
end do
return
end
subroutine gbound(g,tw,w,n,m)
real g(0:8,0:n,0:m)
real w(0:8),tw
! Boundary conditions
! West boundary condition, T=1.
do j=0,m
g(1,0,j)=tw(1-(j/m))(w(1)+w(3))-g(3,0,j)
g(5,0,j)=tw(1-(j/m))(w(5)+w(7))-g(7,0,j)
g(8,0,j)=tw(1-(j/m))(w(8)+w(6))-g(6,0,j)
!g(1,0,j)=tw((j/m))(w(1)+w(3))-g(3,0,j)
!g(5,0,j)=tw((j/m))(w(5)+w(7))-g(7,0,j)
!g(8,0,j)=tw((j/m))*(w(8)+w(6))-g(6,0,j)
! g(1,0,j)=tw*(w(1)+w(3))-g(3,0,j)
!g(5,0,j)=tw*(w(5)+w(7))-g(7,0,j)
!g(8,0,j)=tw*(w(8)+w(6))-g(6,0,j)
end do
! East boundary condition, T=0.
do j=0,m
g(6,n,j)=-g(8,n,j)
g(3,n,j)=-g(1,n,j)
g(7,n,j)=-g(5,n,j)
end do
! Top boundary conditions, Adiabatic
do i=0,n
!g(8,i,m)=g(8,i,m-1)
!g(7,i,m)=g(7,i,m-1)
!g(6,i,m)=g(6,i,m-1)
!g(5,i,m)=g(5,i,m-1)
!g(4,i,m)=g(4,i,m-1)
!g(3,i,m)=g(3,i,m-1)
!g(2,i,m)=g(2,i,m-1)
!g(1,i,m)=g(1,i,m-1)
!g(0,i,m)=g(0,i,m-1)
g(4,i,m)=g(2,i,m)
g(8,i,m)=g(6,i,m)
g(7,i,m)=g(5,i,m)
end do
!Bottom boundary conditions, Adiabatic
do i=0,n
!g(1,i,0)=g(1,i,1)
!g(2,i,0)=g(2,i,1)
!g(3,i,0)=g(3,i,1)
!g(4,i,0)=g(4,i,1)
!g(5,i,0)=g(5,i,1)
!g(6,i,0)=g(6,i,1)
!g(7,i,0)=g(7,i,1)
!g(8,i,0)=g(8,i,1)
!g(0,i,0)=g(0,i,1)
g(2,i,0)=tw*(w(2)+w(4))-g(4,i,0)
g(5,i,0)=tw*(w(5)+w(7))-g(7,i,0)
g(6,i,0)=tw*(w(8)+w(6))-g(8,i,0)
end do
return
end
subroutine tcalcu(g,th,n,m)
real g(0:8,0:n,0:m),th(0:n,0:m)
do j=0,m
do i=0,n
ssumt=0.0
do k=0,8
ssumt=ssumt+g(k,i,j)
end do
th(i,j)=ssumt
end do
end do
return
end
subroutine rhouv(f,rho,u,v,cx,cy,n,m)
real f(0:8,0:n,0:m),rho(0:n,0:m),u(0:n,0:m),v(0:n,0:m),cx(0:8),cy(0:8)
do j=0,m
do i=0,n
ssum=0.0
do k=0,8
ssum=ssum+f(k,i,j)
end do
rho(i,j)=ssum
end do
end do
DO i=0,n
DO j=0,m
usum=0.0
vsum=0.0
DO k=0,8
usum=usum+f(k,i,j)cx(k)
vsum=vsum+f(k,i,j)cy(k)
END DO
u(i,j)=usum/rho(i,j)
v(i,j)=vsum/rho(i,j)
END DO
END DO
return
end
subroutine result(u,v,rho,th,uo,n,m,ra)
real u(0:n,0:m),v(0:n,0:m),th(0:n,0:m)
real strf(0:n,0:m),rho(0:n,0:m)
open(5,file=‘streamt’)
open(6,file=‘nuav’)
! streamfunction calculations
strf(0,0)=0.
do i=0,n
rhoav=0.5(rho(i-1,0)+rho(i,0))
if(i.ne.0) strf(i,0)=strf(i-1,0)-rhoav0.5*(v(i-1,0)+v(i,0))
do j=1,m
rhom=0.5*(rho(i,j)+rho(i,j-1))
strf(i,j)=strf(i,j-1)+rhom0.5(u(i,j-1)+u(i,j))
end do
end do
! ———————————–
write(2,)“VARIABLES =X, Y, U, V”
write(2,)“ZONE “,“I=”,n+1,“J=”,m+1,”,”,“F=BLOCK”
do j=0,m
write(2,)(i/float(m),i=0,n)
end do
do j=0,m
write(2,)(j/float(m),i=0,n)
end do
do j=0,m
write(2,)(u(i,j),i=0,n)
end do
do j=0,m
write(2,)(v(i,j),i=0,n)
end do
do j=0,m
!write(3,)j/float(m),u(5,j)/uo,u(n/2,j)/uo,u(n-10,j)/uo
end do
do i=0,n
write(4,) i/float(n),th(i,m/2)
end do
write(5,)“VARIABLES =X, Y, S, T”
write(5,)“ZONE “,“I=”,n+1,“J=”,m+1,”,”,“F=BLOCK”
do j=0,m
write(5,)(i/float(m),i=0,n)
end do
do j=0,m
write(5,)(j/float(m),i=0,n)
end do
do j=0,m
write(5,)(strf(i,j),i=0,n)
end do
do j=0,m
write(5,)(th(i,j),i=0,n)
end do
! Nusselt number Calculation
snul=0.0
snur=0.0
do j=0,m
rnul=(th(0,j)-th(1,j))float(n)
rnur=(th(n-1,j)-th(n,j))float(n)
snul=snul+rnul
snur=snur+rnur
write(5,)j/float(m),rnul,rnur
end do
avnl=snul/float(m)
avnr=snur/float(m)
write(6,)ra,avnl,avnr
OPEN(21,FILE=‘DATx.m’)
OPEN(22,FILE=‘the.m’)
OPEN(23,FILE=‘DATu.m’)
OPEN(24,FILE=‘DATv.m’)
OPEN(25,FILE=‘steam.m’)
REWIND(21)
REWIND(22)
REWIND(23)
REWIND(24)
REWIND(25)
WRITE(21,)‘x = […’
WRITE(22,)‘the = […’
WRITE(23,)‘u = […’
WRITE(24,)‘v = […’
WRITE(25,)‘stream = […’
DO j=0,m
WRITE(21,)(i,i=0,n)
write(22,100)(th(i,j),i=0,n)
WRITE(23,100)(u(i,j),i=0,n)
WRITE(24,100)(v(i,j),i=0,n)
WRITE(25,100)(strf(i,j),i=0,n)
ENDDO
WRITE(21,*)'] ;'
WRITE(22,*)'] ;'
WRITE(23,*)'] ;'
WRITE(24,*)'] ;'
WRITE(25,*)'] ;'
100 FORMAT(1X,256(E12.6,2X))
return
end]
please correct it if you can