格子玻尔兹曼-恒温无限大板中热扩散D1Q2模型MATLAB计算程序

假设板的初始温度T=O, 当时间t;=0 时，板的左侧面受到高温T=1.0, 板长为100 。
设扩散系数a=0.25, 计算t=200 时，板中的温度分布

1、初始条件：

L=100;%计算区域长度
TL=1;    %温度边界条件
m=100;  %格子数
alpha=0.25; %扩散系数
t=200;  %总时间
dt=L/m;   %时间步长
dx=L/m;   %格子步长
omiga=dt/(alpha+0.5*dt);
f=zeros(t+1,2*(m+1));

2、碰撞步骤

    for j=1:m+1
        T(i,j)=f(i,2*j-1)+f(i,2*j);
        feq(i,2*j-1)=0.5*T(i,j);
        feq(i,2*j)=0.5*T(i,j);
        f(i,2*j-1)=f(i,2*j-1)*(1-omiga)+omiga*feq(i,2*j-1);
        f(i,2*j)=f(i,2*j)*(1-omiga)+omiga*feq(i,2*j);
    end

3、迁移步骤

 for j=1:m
        f(i+1,2*(m+1)-1-2*(j-1))=f(i,2*(m+1)-1-2*j);   
        f(i+1,2*j)=f(i,2*(j+1));  
    end

4、边界条件

    f(i+1,1)=TL-f(i+1,2);
    f(i+1,2*(m+1)-1)=f(i+1,2*(m+1)-3);
    f(i+1,2*(m+1))=f(i+1,2*(m+1)-2); 

5、主循环及绘图

L=100;%计算区域长度
TL=1;    %温度边界条件
m=100;  %格子数
alpha=0.25; %扩散系数
t=200;  %总时间
dt=L/m;   %时间步长
dx=L/m;   %格子步长
omiga=dt/(alpha+0.5*dt);
f=zeros(t+1,2*(m+1));
%主循环
for i=1:t/dt
%       i=4;
    %碰撞步骤
   
    for j=1:m+1
        T(i,j)=f(i,2*j-1)+f(i,2*j);
        feq(i,2*j-1)=0.5*T(i,j);
        feq(i,2*j)=0.5*T(i,j);
        f(i,2*j-1)=f(i,2*j-1)*(1-omiga)+omiga*feq(i,2*j-1);
        f(i,2*j)=f(i,2*j)*(1-omiga)+omiga*feq(i,2*j);
    end

    %迁移步骤
    for j=1:m
        f(i+1,2*(m+1)-1-2*(j-1))=f(i,2*(m+1)-1-2*j);   
        f(i+1,2*j)=f(i,2*(j+1));  
    end

    %边界条件 
    f(i+1,1)=TL-f(i+1,2);
    f(i+1,2*(m+1)-1)=f(i+1,2*(m+1)-3);
    f(i+1,2*(m+1))=f(i+1,2*(m+1)-2); 
    x=[0:dx:L];
    plot(x,T(i,:),'-o')
    xlim([0,30]);
    drawnow;
end