Define velocity vectors In MATLAB can write ex0 10 ex1 11 etc 0 1 2 3 4 5 6 7 8 D2Q9 e 1 e 2 e 3 e 4 e 5 e 6 e 7 e 8 ID: 634834
Download Presentation The PPT/PDF document "Building an LBM Model The following is a..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
Slide1
Building an LBM Model
The following is a mixture of MATLAB and C. Care must be taken, because array indexing begins at 1 in MATLAB and at 0 in C.Slide2
Define velocity vectors
In MATLAB, can write ex(0+1)=0, ex(1+1)=1, etc.
0
1
2
3
4
5
6
7
8
D2Q9
e
1
e
2
e
3
e
4
e
5
e
6
e
7
e
8
%define
e
s
ex(0)= 0;
ey
(0)= 0
ex(1)= 1;
ey
(1)= 0
ex(2)= 0;
ey
(2)= 1
ex(3)=-1;
ey
(3)= 0
ex(4)= 0;
ey
(4)=-1
ex(5)= 1;
ey
(5)= 1
ex(6)=-1;
ey
(6)= 1
ex(7)=-1;
ey
(7)=-1
ex(8)= 1;
ey
(8)=-1Slide3
Problem definition
LY=10
LX=20
tau = 1
g=0.00001
%set solid nodesis_solid_node=zeros(LY,LX)for i=1:LX is_solid_node(1,i)=1 is_solid_node(LY,i)=1end
LX
LY
% if ~is_interior_solid_node(j,i)Slide4
Initialize density and fs (assuming zero velocity)
%define initial density and fs
rho=ones(LY,LX);
f(:,:,1) = (4./9. ).*rho;
f(:,:,2) = (1./9. ).*rho;
f(:,:,3) = (1./9. ).*rho;f(:,:,4) = (1./9. ).*rho;f(:,:,5) = (1./9. ).*rho;f(:,:,6) = (1./36.).*rho;f(:,:,7) = (1./36.).*rho;
f(:,:,8) = (1./36.).*rho;f(:,:,9) = (1./36.).*rho;Slide5
// Computing macroscopic density, rho, and velocity, u=(ux,uy).
for( j=0; j<LY; j++)
{
for( i=0; i<LX; i++)
{ u_x[j][i] = 0.0; u_y[j][i] = 0.0; rho[j][i] = 0.0;
if( !is_solid_node[j][i]) { for( a=0; a<9; a++) { rho[j][i] += f[j][i][a]; u_x[j][i] += ex[a]*f[j][i][a];
u_y[j][i] += ey[a]*f[j][i][a]; } u_x[j][i] /= rho[j][i]; u_y[j][i] /= rho[j][i]; } } }Slide6Slide7
Periodic Boundaries
On boundary, ‘neighboring’ point is on opposite boundary
ip = ( i<LX-1)?( i+1):( 0 );
in = ( i>0 )?( i-1):( LX-1);
jp = ( j<LY-1)?( j+1):( 0 );
jn = ( j>0 )?( j-1):( LY-1);
where
LHS = (COND)?(TRUE_RHS):(FALSE_RHS);
means
if( COND) { LHS=TRUE_RHS;} else{ LHS=FALSE_RHS;}Slide8
// Compute the equilibrium distribution function, feq.
f1=3.;
f2=9./2.;
f3=3./2.; for( j=0; j<LY; j++)
{ for( i=0; i<LX; i++) { if( !is_solid_node[j][i])
{ rt0 = (4./9. )*rhoij; rt1 = (1./9. )*rhoij; rt2 = (1./36.)*rhoij; ueqxij = uxij; //Can add forcing here; see MATLAB code ueqyij = uyij;
uxsq = ueqxij * ueqxij; uysq = ueqyij * ueqyij; uxuy5 = ueqxij + ueqyij; uxuy6 = -ueqxij + ueqyij; uxuy7 = -ueqxij + -ueqyij; uxuy8 = ueqxij + -ueqyij;
usq = uxsq + uysq; Slide9
feqij[0] = rt0*( 1. - f3*usq);
feqij[1] = rt1*( 1. + f1*ueqxij + f2*uxsq - f3*usq);
feqij[2] = rt1*( 1. + f1*ueqyij + f2*uysq - f3*usq);
feqij[3] = rt1*( 1. - f1*ueqxij + f2*uxsq - f3*usq);
feqij[4] = rt1*( 1. - f1*ueqyij + f2*uysq - f3*usq);
feqij[5] = rt2*( 1. + f1*uxuy5 + f2*uxuy5*uxuy5 - f3*usq); feqij[6] = rt2*( 1. + f1*uxuy6 + f2*uxuy6*uxuy6 - f3*usq); feqij[7] = rt2*( 1. + f1*uxuy7 + f2*uxuy7*uxuy7 - f3*usq); feqij[8] = rt2*( 1. + f1*uxuy8 + f2*uxuy8*uxuy8 - f3*usq); }
} }Slide10
// Collision step.
for( j=0; j<LY; j++)
for( i=0; i<LX; i++)
if( !is_solid_node[j][i])
for( a=0; a<9; a++) { fij[a] = fij[a] +
( fij[a] - feqij[a])/tau; }
CollideSlide11
Bounceback Boundaries
Bounceback boundaries are particularly simple
Played a major role in making LBM popular among modelers interested in simulating fluids in domains characterized by complex geometries such as those found in porous media
Their beauty is that one simply needs to designate a particular node as a solid obstacle and no special programming treatment is required Slide12
Two type of solids:Slide13
Bounceback BoundariesSlide14
Collision (solid nodes): Bounceback
// Standard bounceback.
temp = fij[1]; fij[1] = fij[3]; fij[3] = temp;
temp = fij[2]; fij[2] = fij[4]; fij[4] = temp;
temp = fij[5]; fij[5] = fij[7]; fij[7] = temp;
temp = fij[6]; fij[6] = fij[8]; fij[8] = temp;Slide15
// Streaming step.
for( j=0; j<LY; j++)
{
jn = (j>0 )?(j-1):(LY-1) jp = (j<LY-1)?(j+1):(0 )
for( i=0; i<LX; i++) { if( !is_interior_solid_node[j][i])
{ in = (i>0 )?(i-1):(LX-1) ip = (i<LX-1)?(i+1):(0 ) ftemp[j ][i ][0] = f[j][i][0]; ftemp[j ][ip][1] = f[j][i][1];
ftemp[jp][i ][2] = f[j][i][2]; ftemp[j ][in][3] = f[j][i][3]; ftemp[jn][i ][4] = f[j][i][4]; ftemp[jp][ip][5] = f[j][i][5]; ftemp[jp][in][6] = f[j][i][6]; ftemp[jn][in][7] = f[j][i][7];
ftemp[jn][ip][8] = f[j][i][8]; } } }Slide16
Incorporating Forcing (gravity)
Gravitational acceleration is incorporated in a velocity term
ueqxij = u_x(j,i)+tau*g; %add forcing