data.m

晶格波尔兹曼方法

%无量纲形式，△t=1，△x=1，c=△x/△t=1
tic
clear Channel2D
%管道参数初始化
Len_Channel_2D=36;
Channel_2D_half_Width=8;
Width=Channel_2D_half_Width*2;
Channel2D=ones(Len_Channel_2D,Width);
WN=size(Channel2D,2);
%流体属性
cs2=1/3;
cP_visco=0.5;%动力粘度,单位:cP
density=1.0;
Lky_visco=cP_visco/density;%格子运动粘度 
omega=(Lky_visco/cs2+0.5).^-1;%松弛频率

%初始速度
uy0=0; 
ux0=0;%这个0.0001是怎么定的？

%解析速度（Poiseuille）
dPdL=-0.0125;%压力梯度
uy_fin_max=dPdL*(Channel_2D_half_Width.^2)/(2*Lky_visco);
x_profile=(-Channel_2D_half_Width:Channel_2D_half_Width-1)+0.5;
uy_analy_profile=uy_fin_max.*(1-(x_profile/Channel_2D_half_Width).^2);
%流场平均速度初始化
av_vel_t=1.e+10; 
%管道层数扩展
inb=2;%入口扩展层数
oub=2;%出口扩展层数
inlet=ones(inb,WN);
outlet=ones(oub,WN);
Channel2D=[inlet;Channel2D;outlet];
LN=size(Channel2D,1);
wb=2;%边界扩展层数
Channel2D=[zeros(LN,wb),Channel2D,zeros(LN,wb)];
[LN WN]=size(Channel2D);%最终的整个区域是40*16
%画解析速度曲线
uy_analy_profile=[zeros(1,wb),uy_analy_profile,zeros(1,wb) ];
x_pro_fig=[x_profile(1)-(wb:-1:1),[x_profile,(1:wb)+x_profile(end)]];
figure(20);plot(x_pro_fig,uy_analy_profile,'-');grid on;
title('Analytical parab. profile for Poiseuille planar flow in a channel')
%计算内部，边界和整体区域（40*16）索引号
Channel2D=logical(Channel2D);
Obstacles=bwperim(Channel2D,8);
border=true([LN,WN]);
border([1:inb,LN-oub:LN],wb+2:WN-wb-1)=0;
Obstacles=Obstacles.*border;
[iWhole jWhole]=find(Channel2D==1);lenWhole=length(iWhole);ijWhole=(jWhole-1)*LN+iWhole;%整体区域
[iObs jObs]=find(Obstacles);lenObs=length(iObs);ijObs=(jObs-1)*LN+iObs;%边界区域
[iInner jInner]=find((Channel2D==1&~Obstacles));lenInner=length(iInner);ijInner=(jInner-1)*LN+iInner;%内部区域
%方向信息
%  6 2 5        NW  N  NE
%  3 9 1 ---->  W   RP  E
%  7 4 8        SW  S  SE
East=1;North=2;West=3;South=4;NE=5;NW=6;SW=7;SE=8;RP=9;
C_x=[1 0 -1  0 1 -1 -1  1 0]; 
C_y=[0 1  0 -1 1  1 -1 -1 0]; 
C=[C_x;C_y];
N_c=9;%方向数
ic_op = [3 4 1 2 7 8 5 6];%速度反向
yi2=[LN,1:LN,1];%周期边界条件
%平衡函数系数
w0=4/9;w1=1/9;w2=1/36;
W=[w1 w1 w1 w1 w2 w2 w2 w2 w0];
f1=3.0;f2=4.5; f3=1.5;
%变量定义
NxM=LN*WN;
F=zeros(LN,WN,N_c);
F_EQ=zeros(LN,WN,N_c);
rho=ones(LN,WN);
temp1=zeros(LN,WN);
ux=zeros(LN,WN);uy=zeros(LN,WN);uyout=zeros(LN,WN);%无量纲速度
uxsq=zeros(LN,WN);uysq=zeros(LN,WN);usq=zeros(LN,WN);%速度平方
% 分布函数初始化
for ia=1:lenWhole
	i=iWhole(ia);
	j=jWhole(ia);
	F(i,j,:)=1/9;
end

%速度初始化
uy(ijWhole)=uy0;
ux(ijWhole)=ux0;
rho(ijWhole)=density;
force=dPdL/6*[0 1 0 -1 1 1 -1  -1  0]';%[0 -1 0 1 -1 -1 1  1  0] ---> E  N E S NE NW SW SE RP

%控制参数
StopFlag=false;
Max_Iter=3000;
Check_Iter=1; 
Output_Every=20;
Cur_Iter=0;
toler=1.0e-8;
Cond_path=[];
density_path=[];