d2q9.m

晶格波尔兹曼方法

Data;
while(~StopFlag)
    Cur_Iter=Cur_Iter+1;
    %宏观量计算
    rho=sum(F,3);                         %我的理解：得到每个格点的局部密度
    
    if Cur_Iter >1                        %有什么必要么？
        ux=zeros(LN,WN);
        uy=zeros(LN,WN);
        for ic=1:N_c-1;                   %求出了宏观动量
            ux = ux + C_x(ic).*F(:,:,ic);
            uy = uy + C_y(ic).*F(:,:,ic);
        end
    end
    ux(ijWhole)=ux(ijWhole)./rho(ijWhole);
    uy(ijWhole)=uy(ijWhole)./rho(ijWhole);%求出了宏观速度
    
    uxsq(ijWhole)=ux(ijWhole).^2;
    uysq(ijWhole)=uy(ijWhole).^2; 
    usq(ijWhole)=uxsq(ijWhole)+uysq(ijWhole);%求出了个格点宏观速度的平方，为计算平衡分布函数做准备
    %***************************************
    %          粒子碰撞（整个区域）
    %***************************************
    %方向信息
    %  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;
    rt0=w0.*rho;rt1=w1.*rho;rt2=w2.*rho;
    F_EQ(ijWhole+NxM*(1-1))=rt1(ijWhole).*(1+f1*ux(ijWhole)+f2*uxsq(ijWhole)-f3*usq(ijWhole));
    F_EQ(ijWhole+NxM*(9-1))=rt0(ijWhole).*(1-f3*usq(ijWhole));
    F_EQ(ijWhole+NxM*(2-1))=rt1(ijWhole).*(1+f1*uy(ijWhole)+f2*uysq(ijWhole)-f3*usq(ijWhole));
    F_EQ(ijWhole+NxM*(3-1))=rt1(ijWhole).*(1-f1*ux(ijWhole)+f2*uxsq(ijWhole)-f3*usq(ijWhole));
    F_EQ(ijWhole+NxM*(4-1))=rt1(ijWhole).*(1-f1*uy(ijWhole)+f2*uysq(ijWhole)-f3*usq(ijWhole));
    F_EQ(ijWhole+NxM*(5-1))=rt2(ijWhole).*(1+f1*(+ux(ijWhole)+uy(ijWhole))+f2*(+ux(ijWhole)+uy(ijWhole)).^2-f3.*usq(ijWhole));
    F_EQ(ijWhole+NxM*(6-1))=rt2(ijWhole).*(1+f1*(-ux(ijWhole)+uy(ijWhole))+f2*(-ux(ijWhole)+uy(ijWhole)).^2-f3.*usq(ijWhole));
    F_EQ(ijWhole+NxM*(7-1))=rt2(ijWhole).*(1+f1*(-ux(ijWhole)-uy(ijWhole))+f2*(-ux(ijWhole)-uy(ijWhole)).^2-f3.*usq(ijWhole));
    F_EQ(ijWhole+NxM*(8-1))=rt2(ijWhole).*(1+f1*(+ux(ijWhole)-uy(ijWhole))+f2*(+ux(ijWhole)-uy(ijWhole)).^2-f3.*usq(ijWhole)); 
    F=(1.-omega).*F+omega.*F_EQ;%注意omega的具体意义
    %在流体上施加压力
    for ic=1:N_c                              %外层为ic方向循环，内层为格子点循环
        for ia=1:lenWhole
            i=iWhole(ia);  
            j=jWhole(ia);
            F(i,j,ic)= F(i,j,ic) + force(ic); %force已经在data.m中定义好了，注意这是常力矩阵，我的需要在这个地方做文章！
        end
    end
    %***************************************
    %               粒子迁移
    %***************************************   
    F_EQ = F;
    for ic=1:1:N_c-1       %碰撞前的速度向量
        ic2=ic_op(ic);     %碰撞后的速度向量，利用前一个当索引
        temp1=F_EQ(:,:,ic);%存储了每个方向上所有点的平衡分布函数
        %内部区域
        for ia=1:1:lenInner
            i=iInner(ia);
            j=jInner(ia);
            i2=i+C_y(ic);        %由于det_x=1,det_t=1,所以i,j既是索引，又可以代表空间的位置坐标
            j2=j+C_x(ic);
            i2=yi2(i2+1);        %周期性边界条件，巧妙，服了！但是将inlet,outlet设成两层的意义是什么？
            F(i2,j2,ic)=temp1(i,j);%将该方向上的分布函数迁移到运动后的位置
        end
        %边界                    %只有边界层的格子才有必要验证是否用反弹边界条件
        for ia=1:1:lenObs
            i=iObs(ia);
            j=jObs(ia);
            i2=i+C_y(ic);
            j2=j+C_x(ic);        
            i2=yi2(i2+1);        %竖直方向仍然要考虑周期条件
            if Channel2D(i2,j2)==0
                F(i,j,ic2) =temp1(i,j);%如果碰到管壁反弹回去，反向增加分布函数
            else
                F(i2,j2,ic)=temp1(i,j);
            end
        end
    end 
    
    
    %重新计算宏观速度和密度
    rho=sum(F,3);
    uyout= zeros(LN,WN);
    for ic=1:N_c-1;
        uyout= uyout + C_y(ic).*F(:,:,ic);
    end
    %收敛条件计算
    if mod(Cur_Iter,Check_Iter)==0
        vect=rho(ijWhole); 
        cur_density=mean(vect);
        vect=uy(ijWhole);
        av_vel_int=mean(vect);
        av_vel_tp1=av_vel_int; 
        Condition=abs(abs(av_vel_t/av_vel_tp1 )-1);%收敛系数
        Cond_path=[Cond_path,Condition];%记录收敛系数
        density_path=[density_path,cur_density];%记录密度
        av_vel_t=av_vel_tp1;
        if (Condition < toler)||(Cur_Iter > Max_Iter)
            StopFlag=true;
            display('Stop iteration: Convergence met or iteration exeeding the max allowed')
            display(['Current iteration: ',num2str(Cur_Iter),' Max Number of iter: ',num2str(Max_Iter)])
            break;
        end
    end
    %输出
    if mod(Cur_Iter,Output_Every)==0
        rho=sum(F,3);
        up=2; %从较低层可视化
        figure(15);hold off;feather(ux(LN-up,:),uy(LN-up,:));
        figure(15);hold on;plot(uy_analy_profile,'r-');
        title('Analytical vs LB calculated, fluid velocity parabolic profile');
        pause(3);
    end
end 
figure;plot(Cond_path(2:end));title('convergence path');
figure;plot(uy(LN-up,:)-uy_analy_profile);title('difference : LB - Analytical solution');
toc