感谢群里的朋友们提供帮助！还是老样子，有啥问题Feel free to tell us~毕竟群众力量大嘛~格子玻尔兹曼救星QQ群：293267908。

         流体计算领域中，LBM还是个比较新的思想，最近宝宝正在尝试性的进行研究。首先说一下教材吧。我在某宝买了中文版，杨大勇翻译的，中国工信出版。 还可以看英文原版哈，英文教材链接：https://pan.baidu.com/s/1D6JFR9zp7DGfduvb-4xSUQ，提取码：s5cr。

     开博客一方面自己做个学习记录，也算方便大家看书了。因为第一次用F所以加了很多自己的备注，大家看得难受忽略就好啊。不足之处欢迎大家一同学习。

---------------------------------------------------------------------

     第一章巴拉巴拉的其实说了以下3个事：

     1.1 Boltzmann方程是介观尺度。

     1.2 宏观表现出的温度与压力正比的关系，物理本质上是由于粒子动能与温度正比。

     1.3 所谓分布函数f(c)，指的是多个粒子在某速度范围内（dc）的概率。当流体处于热平衡状态时，粒子动量的分布函数是不变的。

     第二章 Boltzmann方程又是啥呢？

     2.1 分布函数变成平衡分布函数之前，总是会变，变的原因是碰撞。分布函数的总变化率就是速率。

     2.2 既然要求解分布函数，就不得不说怎么求碰撞。用平衡分布函数-分布函数，就是变化量的净值，然后除以松弛因子。我个人理解松弛因子是用以控制到达平衡点的速度，即每计算一个时间周期，系统由不平衡状态向平衡状态前进（feq-f）/τ。

     2.3 D1Q3就是1威三个速度矢量，以此类推。

    第三章 扩散

     3.3 分布函数的总和和平衡分布函数的总和都等于因变量。即，流体熵增加，高动能的粒子动能减少而低动能的粒子动能增加，但是总值不变。所以无论是否达到最终平衡，所有动能的总和是不变的，所以可以看到分布函数的总和和平衡分布函数的总和都是因变量。但是feq和f不一定相等。一旦相等了，就是平衡状态了。

     3.5  展开的目的，是获得宏观和介观尺度的关系。

      3.5.2中所谓的不保存前值f(x,t)和为了避免覆盖新值，是因为如果每个点都同时向左向右迁移，会相互覆盖。

   第五章 举个栗子：方腔流。有错误的地方我做了标记。

       call的使用还是要注意输入顺序的。比如collesion(f,feq,rho,w,cy,cx,u,v,omega,n,m)，括号里面的引用变量的顺序，子函数和主函数中必须保持顺序一致。

   ! D2Q9,书P141.   等温不可压缩，lid-driven cavity
   !英文版：     《Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes》
  
   
   !设置参数
    parameter (n=100,m=100)          !设置格子数
    real f(0:8,0:n,0:m)                      !分布函数
    real feq(0:8,0:n,0:m),rho(0:n,0:m)  !局部平衡分布函数
    real w(0:8),cx(0:8) ,cy(0:8)             !0~8，共九个方向的权重因子w(0:8)；沿着格子迁移方向的单位矢量c，分别具有xy方向的单位矢量的分量。
    real u(0:n,0:m),v(0:n,0:m)              !xy方向速度
    integer i                                      !英文版有错。???????
    open(2,file='velocity100points.txt')                   !.xls也可以
    open(3,file='uvely.txt')
    open(4,file='vvelx.txt')
    open(8,file='timeu.txt')
    !----------
    uo=0.10                                                !初始速度       
    sumvelo=0.00                                   !?
    rhoo=5.00                                              !密度
    dx=1.0
    dy=dx
    dt=1.0
    alpha=0.01                                             !黏度，在传热中是扩散系数
    Re=uo*m/alpha                                       !雷诺数
    print *,"Re=",Re
    ck=dx/dt
    csq=ck*ck                                               !注意，此处不同于书上。但是当dt=dx时，csq写了和没写计算值一样的
    omega=1.0/(3.*alpha/(csq*dt)+0.5)         !中文书71页
    mstep=400                                             !步数，可以先设置小一点！
    w(0)=4./9.                                              !抄书就行
do i=1,4
w(i)=1./9.
end do
do i=5,8
w(i)=1./36.
end do
cx(0)=0
cx(1)=1
cx(2)=0
cx(3)=-1
cx(4)=0
cx(5)=1
cx(6)=-1
cx(7)=-1
cx(8)=1
cy(0)=0
cy(1)=0
cy(2)=1
cy(3)=0
cy(4)=-1
cy(5)=1
cy(6)=1
cy(7)=-1
cy(8)=-1
        do j=0,m
        do i=0,n
            rho(i,j)=rhoo                                 !因变量。热扩散中是温度，动量扩散中是速度，质量扩散中是组分数。在这里初速度为0，所以因变量应该是密度吧？？
             u(i,j)=0.0
             v(i,j)=0.0
        end do
        end do
        do i=1,n-1      !为了配合原著，此处暂时修改为n-1                                            !原书此处i只有n-1个值！实际是一百个节点。改为n。
            u(i,m)=uo                                      !设置上端速度的初始值，uo为x轴正方向=0.10
            v(i,m)=0.0
        end do
    !主程序开始
1 do kk=1,mstep                                 !原书此处有误。已改
        call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
        call streaming(f,n,m)
        call sfbound(f,n,m,uo)
        call rhouv(f,rho,u,v,cx,cy,n,m)
            write(8,*)  kk,u(n/2,m/2),v(n/2,m/2)            !timeu!!!!!!
            !print *,"i=",i, "j=",j,",","u(i,j)=",u(i,j),",","v(i,j)=",v(i,j),",","rho(i,j)=",rho(i,j)
                          !发现问题：ij最大值超过100；速度均为0！！
        end do
    !主循环完
        call result(u,v,rho,uo,n,m)
        stop
        end
    !主程序完

subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)   !计算平衡分布指数
   real f(0:8,0:n,0:m)                                                 !计算平衡分布指数之前，把参数范围写一下
   real feq(0:8,0:n,0:m),rho(0:n,0:m)
   real w(0:8), cx(0:8),cy(0:8)
   real u(0:n,0:m), v(0:n,0:m)
   DO i=0,n
   DO j=0,m
   t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
        do k=0,8
            t2=u(i,j)*cx(k)+v(i,j)*cy(k)                               !第一次编辑，我在这里错了！已改正！
            feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)   !中文书70页。平衡函数
            f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)       !扩散过程的每一步都要乘以松弛频率omega
                     !print *,"k=",k,"i=",i, "j=",j,"f(k,i,j)=",f(k,i,j),"feq(k,i,j)=",feq(k,i,j)     !此处用于检测
        end do
        end do
        end do
 return
 end

subroutine streaming(f,n,m)                    !碰撞和扩散
    real f(0:8,0:n,0:m)
    ! 13之间,
        do j=0,m
        do i=n,1,-1
            f(1,i,j)=f(1,i-1,j)              !从右到左
                        ! print *,"i=",i,"j=",j,"f(1,i,j)=",f(1,i,j)        !此处用于检测1
        end do
        do i=0,n-1
            f(3,i,j)=f(3,i+1,j)             !从左到右
                        ! print *,"i=",i,"j=",j,"f(3,i,j)=",f(3,i,j)        !此处用于检测2
        end do
        end do
    !!!!!!!!!!!!
        do j=m,1,-1                         !从上到下
        do i=0,n
            f(2,i,j)=f(2,i,j-1)
                        ! print *,"i=",i,"j=",j,"f(2,i,j)=",f(2,i,j)        !此处用于检测3
        end do
        do i=n,1,-1
            f(5,i,j)=f(5,i-1,j-1)     
                        !  print *,"i=",i,"j=",j,"f(5,i,j)=",f(5,i,j)        !此处用于检测  4   
        end do
        do i=0,n-1
            f(6,i,j)=f(6,i+1,j-1)      
                         ! print *,"i=",i,"j=",j,"f(6,i,j)=",f(6,i,j)        !此处用于检测 5（i=           1 j=         100 f(6,i,j)=  0.2620545    
        end do
        end do
    !!!!!!!!!!!!
        do j=0,m-1                         !从下到上
        do i=0,n
            f(4,i,j)=f(4,i,j+1)
                         ! print *,"i=",i,"j=",j,"f(4,i,j)=",f(4,i,j)        !此处用于检测 6     
        end do
        do i=0,n-1
            f(7,i,j)=f(7,i+1,j+1)          
                         ! print *,"i=",i,"j=",j,"f(7,i,j)=",f(7,i,j)        !此处用于检测 7      
        end do
        do i=n,1,-1
            f(8,i,j)=f(8,i-1,j+1)     
                         !  print *,"i=",i,"j=",j,"f(8,i,j)=",f(8,i,j)        !此处用于检测 8（i=         100 j=           0 f(8,i,j)=  0.2620545  i=          99 j=           0 f(8,i,j)=  0.2620545）  !以上均正确！         
        end do
        end do
    return
    end

subroutine sfbound(f,n,m,uo)
        real f(0:8,0:n,0:m)
        do j=0,m
        !west
             f(1,0,j)=f(3,0,j)
             f(5,0,j)=f(7,0,j)
             f(8,0,j)=f(6,0,j)
                    !print *,"j=",j,"f(8,0,j)=",f(8,0,j)    !此处用于检测9（ j=           0 f(8,0,j)=  0.2620545）
        !east
             f(3,n,j)=f(1,n,j)
             f(7,n,j)=f(5,n,j)
             f(6,n,j)=f(8,n,j)
                     !print *,"j=",j,"f(6,n,j)=",f(6,n,j)    !此处 用于检测10（
        end do
        do i=0,n
        !south
             f(2,i,0)=f(4,i,0)
             f(5,i,0)=f(7,i,0)
             f(6,i,0)=f(8,i,0)
                   !print *,"i=",i,"f(6,i,0)=",f(6,i,0)       !此处用于检测11
        end do
        !north   注意北侧有运动
             do i=1,n-1
             rhon=f(0,i,m)+f(1,i,m)+f(3,i,m)+2.*(f(2,i,m)+f(6,i,m)+f(5,i,m))
             f(4,i,m)=f(2,i,m)
             f(7,i,m)=f(5,i,m)-rhon*uo/6.0   !此处与79页不符，注意公式 f1=f3+2/3(ρu)!!!!!!!!!!!!!!!!中文书（5.21）我第一次没看到郁闷了两天！
             f(8,i,m)=f(6,i,m)+rhon*uo/6.0   !
                  !print *,"i=",i, "f(8,i,0)=",f(8,i,0)         !此处用于检测12（）
             end do
             return
             end

subroutine rhouv(f,rho,u,v,cx,cy,n,m)    ! 这里有问题！！！
    real f(0:8,0:n,0:m)
    real rho(0:n,0:m)
    real cx(0:8) ,cy(0:8)
    real u(0:n,0:m),v(0:n,0:m)
                        !print *, "cx(3)=",cx(3),"cy(5)=",cy(5)          !此处发现问题！ 13（cx(3)=  -1.000000     cy(5)=   1.000000）
        do j=0,m
        do i=0,n
            ssum=0.0
        do k=0,8
            ssum=ssum+f(k,i,j)
            !print *,"k=",k,"i=",i, "j=",j, "ssum=",ssum               !此处用于检测 14
        end do
            rho(i,j)=ssum
        end do
        end do
    do i=1,n
        rho(i,m)=f(0,i,m)+f(1,i,m)+f(3,i,m)+2.*(f(2,i,m)+f(6,i,m)+f(5,i,m))
        !print *,"i=",i,"rho(i,m)=",rho(i,m)                                !此处用于检测！15
    end do
    do i=1,n
    do j=1,m-1                      !留意！是否超百？？？
        usum=0.0
        vsum=0.0
    do k=0,8
        usum=usum+f(k,i,j)*cx(k)                                            !这个方向的动量等于分布函数乘以在这个方向的速度分量
        vsum=vsum+f(k,i,j)*cy(k)                                            !此处第一次书写就写错了！
                     !print *, "k=",k,"i=",i, "j=",j
                     !print *,"f(k,i,j)=",f(k,i,j)                             !此处用于检测 16（）
                     !print *, "cx(k)=",cx(k),"cy(k)=",cy(k)          !此处用于检测    发现问题c等于0！！ 17（）
    end do
                     !print *, "usum=",usum,"vsum=",vsum        !此处用于检测！   usum=  0.0000000E+00 vsum=  0.0000000E+00
        u(i,j)=usum/rho(i,j)                                                    !此处第一次书写就写错了!
        v(i,j)=vsum/rho(i,j)                                                    !动量除以质量（单位体积的质量其实就是速度），结果就是速度了
!                     print *,"rho(i,j)=",rho(i,j)                         !此处用于检测     18（）
!                     print *,"u(i,j)=",u(i,j)                               !此处用于检测    19（）发现问题！！！u、vsum等于0导致的全局速度为0！！  
!                     print *,"v(i,j)=",v(i,j)                               !此处用于检测    发现问题！！！u、vsum等于0导致的全局速度为0！！
    end do
    end do
    return
    end
    !---------------
    !---------------结果输出
subroutine result(u,v,rho,uo,n,m)
    real u(0:n,0:m),v(0:n,0:m)
    real rho(0:n,0:m),strf(0:n,0:m)    
    !---------------随后的内容为自己写
    !real rho(0:n,0:m),strfv(0:n,0:m)                              !此处出现错误！！！！！！！！！！！！！！
open(5,file='streamf.txt')  !动量部分
    strf(0,0)=0.0
  do i=1,n   
        rhoav=0.5*(rho(i-1,0)+rho(i,0))                                            !此处原书出现x=-1，已修改
        if(i.ne.0) strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))          !?为什么strf-rhoav？
                                                                                                 !为什么在x轴上对y方向的速度进行计算？？
     do j=1,m
                rhom=0.5*(rho(i,j)+rho(i,j-1))
                strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j)+u(i,j-1))                  !为什么只对x方向速度求解？ !?为什么strf+rhom？
     end do
  end do

      !------
write(2,*)"VARIABLES=X,Y,U,V,S"   ! uvfield
write(2,*)"ZONE","，","I=",n,"，","J=",m,"，","F=BLOOK"
do j=0,20  
    do i=0,20
        write(2,*),i,j,u(i,j),v(i,j)
    end do
end do

!    do j=0,m  
!        write(2,*)(i,i=0,n)   由于要使用tecplot，我重新写了这一块，并输出为:uvfield_tecplot.txt
!    end do
!    do j=0,m
!        write(2,*)(j,i=0,n)     
!    end do
!    do j=0,m
!        write(2,*)(u(i,j),i=0,n)
!    end do
!    do j=0,m
!        write(2,*)(v(i,j),i=0,n)
!    end do
!    do j=0,m
!        write(2,*)(strf(i,j),i=0,n)
!    end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    do j=0,m
        write(3,*)j/float(m),u(n/2,j)/uo,u(n/4,j)/uo,u(3*n/4,j)/uo   !uvely
    end do   
    do i=0,n
        write(4,*)i/float(n),v(i,m/2)/uo                                         !vvelx
    end do    
return
end

%流线图MATLAB源程序如下：
close all;clc;clear all;

text= load('velocity10000POINTS.txt');
x=text(:,1);
r=text(:,2);
dvx=text(:,3);
dvr=text(:,4);

%%
Fx = scatteredInterpolant(x,r,dvx);  %对数据集执行插值
Fr = scatteredInterpolant(x,r,dvr);
%%
xx=linspace(min(x),max(x),600);      % xx= linspace(x1,x2,n) 生成 n 个点。这些点的间距为 (x2-x1)/(n-1)。 调节此处可以调整疏密度！
rr=linspace(min(r),max(r),600);
[xgg,rgg]=meshgrid(xx,rr);
xstream = Fx(xgg,rgg);
ystream = Fr(xgg,rgg);
%%
scrsz = get(0,'ScreenSize');         %得到屏幕参数

figure1 = figure('Position',[0.05*scrsz(3) 0.05*scrsz(4) 0.9*scrsz(3) 0.9*scrsz(4)]);   % 改变画图大小位置
                                     % scrsz(1): 屏幕最左坐标；scrsz(2): 屏幕最下坐标
                                     % scrsz(3): 屏幕宽(像素)；% scrsz(4): 屏幕高(像素)
[xs,rs] = meshgrid(x,r);
[dvxs,dvrs] = meshgrid(dvx,dvr);
quiver(x,r,dvx,dvr,'r');

numstream=900;
strx=randi([2,99],numstream,1);       %randi([imin,imax],...) 返回一个在[imin,imax]范围内的伪随机整数
stry=randi([0,40],numstream,1);       %r = randi(imax,[m,n])，返回一个在[1,imax]范围内的的m*n的伪随机整数矩阵 原始为stry=randi([0,12],numstream,1);
                                      %值strx、y代表流线的起始位置
strx=[strx,strx];
stry=[stry,-stry];

h=streamline(xgg,rgg,xstream,ystream,strx,stry); %streamline(X,Y,Z,U,V,W,startx,starty,startz) 根据三维向量数据 U、V 和 W 绘制流线图。
                                                 %数组 X、Y 和 Z 用于定义 U、V 和 W 的坐标，它们必须是单调的，无需间距均匀。X、Y 和 Z 必须具有相同数量的元素，就像由 meshgrid 生成一样。
                                                 %startx、starty 和 startz 定义流线图的起始位置。
set(h,'LineWidth',0.5,'Color','k')

axis equal
axis tight
box on

%% 2nd
figure2 = figure('Position',[0.05*scrsz(3) 0.05*scrsz(4) 0.9*scrsz(3) 0.9*scrsz(4)]);

[xs,rs] = meshgrid(x,r);
[dvxs,dvrs] = meshgrid(dvx,dvr);
quiver(x,r,dvx,dvr,'r');

numstream=800;
strx=randi([2,99],numstream,1);
stry=randi([0,40],numstream,1);         %此处原为[0,12]

strx=[strx,strx];
stry=[stry,-stry];

h=streamslice(xgg,rgg,xstream,ystream);  %流线图

set(h,'LineWidth',0.7,'Color','b')

axis equal
axis tight
box on

科研党有问题可以联系QQ邮箱[email protected]一同进步.或加QQ群293267908。