%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