LBM方法

LBM方法，多弛豫时间的重力驱动poiseuille流的matlab代码

%GLBE based incompressible MRT model for gravity 
%driven fluid flow simulation in a parallel channel.
%
%D2Q9: 2D- 9 velocity model

%lx is length of domain and ly is width. 
%g is applied gravity in X-direction
%Copyright: Shadab Anwar Florida International University USA.
%sanwa001@fiu.edu
%*********************************************************************
clear all;close all;
wa（9）=4/9;wa（1:4）=1/9;wa（5:8）=1/36;
lx=8;ly=33; rho0_in=1.0001; rho0_out=1.;
u=zeros（lylx2）; g= 1e-04;
u_pois=zeros（1ly）; f=zeros（lylx9）;
ftemp=zeros（lylx9）;rho=ones（lylx）;
u（::1）=0.0041;
M=[   111111111;
     -1-1-1-12222-4;
     -2-2-2-211114;
      10-101-1-110;
     -20201-1-110;
      010-111-1-10;
      0-20211-1-10;
      1-11-100000;
      00001-11-10;
];

s=[0-1.64-1.540-1.90-1.9-1.-1.]; % relaxation parameters
%s（9） controls the kinematic viscosity and hence the Reynolds number.

ts=300;

a2=-8;a3=4;c1=-2;g1=2/3;g3=2/3;g2=18;g4=-18;  % a3 and g4 are free to tune
e（1:）=[10-101-1-110]; %ex
e（2:）=[010-111-1-10]; %ey

%i=10;j=5;
%rho=1;

 tau（1）=-1/s（9）; 
 nu=（1/3）*（tau（1）-0.5）;
 ni=ly;
 %dpdl=（1/3）*（1e-04）/lx;
  mm=1;

for i=-（ni-1）/2:（ni-2）/2
u_pois（mm）=（g/（2*nu））*（（（ni-2）/2）^2-i^2）;
%u_pois（mm）=（dpdl/（2*nu））*（（（ni-2）/2）^2-i^2）;
u_pois（1）=0.0;
mm=mm+1;
end
u_pois（mm）=0.0;

%initialize f‘s
for i=1:lx
    for j=1:ly

u（ji1）=u_pois（j）;

jx=u（ji1）;
jy=u（ji2）;

m_eq（1） = rho（ji）; %density
m_eq（2） = （a2/4）*rho（ji） + （1/6）*g2*（jx^2+jy^2）;  %energy
m_eq（3） = （a3/4）*rho（ji） + （1/6）*g4*（jx^2+jy^2）;  %energy square
m_eq（4） = jx; %momentum in x-dir
m_eq（5） = （1/2）*c1*jx; %energy flux in x-dir
m_eq（6） = jy; %momentum in y-dir
m_eq（7） = （1/2）*c1*jy; %energy flux in y-dir
m_eq（8） = （3/2）*g1*（jx^2 - jy^2）; %diagonal comp of stress tensor
m_eq（9） = （3/2）*g3*jx*jy;  %off diagonal comp of stress tensor
            
f（ji:）=inv（M）*（m_eq‘）;
    end
end

        F（9） = 3*wa（9）*（ e（19）*g ）;
        for a=1:8
             F（a）=3*wa（a）*（ e（1a）*g ）;  %F is applied body force （gravity）
        end
for t=1:ts
    t
%Macroscopic variable
for i=1:lx
    for j=1:ly
        
        rho（ji）=sum（f（ji:））;
        u（ji:）=0.;
        
         for a=1:9
        u（ji1） = u（ji1） + e（1a）*f（jia）;
        u（ji2） = u（ji2） + e（2a）*f（jia）;
          end
          
      u（ji1）= u（ji1）/rho（ji）;
      u（ji2）= u（ji2）/rho（ji）;
  
f_pre=ones（91）;f_post=ones（91）;

%Compute Meq or moment of feq
jx=u（ji1）;
jy=u（ji2）;

m_eq（1） = rho（ji）; %density
m_eq（2） = （a2/4）*rho（ji） + （1/6）*g2*（jx^2+jy^2）;  %energy
m_eq（3） = （a3/4）*rho（ji） + （1/6）*g4*（jx^2+jy^2）;  %energy square
m_eq（4） = rho（ji）*jx; %momentum in x-dir
m_eq（5） = （1/2）*c1*jx; %energy flux in x-dir
m_eq（6） = rho（ji）*jy; %momentum in y-dir
m_eq（7） = （1/2）*c1*jy; %energy flux in y-dir
m_eq（8）

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        4724  2013-04-28 18:45  Untitled.m