格子玻尔兹曼法学习记录（附MATLAB画图源程序）
感谢群里的朋友们提供帮助！还是老样子，有啥问题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
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
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)
end do
do i=n,1,-1
f(5,i,j)=f(5,i-1,j-1)
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) !此处用于检测 （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)
end do
do i=0,n-1
f(7,i,j)=f(7,i+1,j+1)
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) !此处用于检测 （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)
!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)
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)
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)!!!!
f(8,i,m)=f(6,i,m)+rhon*uo/6.0 !
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)
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
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)
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 *, "cx(k)=",cx(k),"cy(k)=",cy(k) !此处用于检测 发现问题c等于0！！ 17（）
end do
u(i,j)=usum/rho(i,j)
v(i,j)=vsum/rho(i,j) !动量除以质量（单位体积的质量其实就是速度），结果就是速度了
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; % number of lines
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); %流线图
%h=streamslice(xgg,rgg,xstream,ystream); %流线图2, with arrows
set(h,'LineWidth',0.7,'Color','b')
axis equal
axis tight
box on
科研党有问题可以联系QQ邮箱1542170984@qq.com一同进步.或加QQ群293267908。
格子玻尔兹曼法学习记录（附MATLAB画图源程序）相关推荐
- 格子玻尔兹曼流体代码_格子玻尔兹曼方法（LBM）学习：对流-扩散问题（附MATLAB代码）...
(๑❛ᴗ❛๑) 麻烦各位读者收藏之余点个喜欢或赞呢,咱也更有干劲了~ OrzSunspot:格子玻尔兹曼方法(LBM)学习:等温不可压缩流体流动问题(附MATLAB代码)zhuanlan.zhihu ...
- 格子玻尔兹曼(LBM)小白的进阶之路
格子玻尔兹曼(LBM)小白的进阶之路 起因 2020年4月16日,距离新型冠状病毒爆发已经五个月的时间了,作为一名科研人员,高校迟迟未发开学通知,作为一名实验人员,第一次那么迫切希望可以熬夜做实验.但 ...
- python可以应用lbm_格子玻尔兹曼方法（LBM）python程序提速
研究生开学已经两周了,一直在学习跟LBM相关的编程知识.由于自己数值传热学的基础不是太好,为了能够快速地融入到现有的工作当中我将工作重心侧重在了编程方面,而不是相关模型和边界条件等的学习.我的主要参考 ...
- 格子玻尔兹曼在多孔介质孔隙尺度气泡输运调控中的应用和MATLAB仿真『需要数据和代码请先私信』
引言 格子玻尔兹曼方法(Lattice Boltzmann Method, LBM)是一种基于介观(mesoscopic)模拟尺度的计算流体力学方法.该方法相比于其他传统CFD计算方法,具有介于微观分 ...
- 格子玻尔兹曼方法书中，计算机代码（Fortran语言）FDM的输出结果是什么，为什么显示程序“[25024] Console1.exe”已退出，返回值为 0 (0x0)。
格子玻尔兹曼方法书中,计算机代码(Fortran语言)FDM的输出结果是什么,为什么显示程序"[25024] Console1.exe"已退出,返回值为 0 (0x0).
- matlab lbm 代码,Matlab实现格子玻尔兹曼方法(Lattice Boltzmann Method,LBM)模拟
%LBM的matlab代码 %Matlab实现格子玻尔兹曼方法(Lattice Boltzmann Method,LBM)模拟 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % c ...
- matlab boltzmann函数,Matlab实现格子玻尔兹曼方法(Lattice Boltzmann Method,LBM)模拟
Matlab实现格子玻尔兹曼方法(Lattice Boltzmann Method,LBM)模拟 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % cylinder.m: Flow ...
- 求助 格子玻尔兹曼数值模拟
需要将一份用于圆柱扰流的换热格子玻尔兹曼程序耦合一个温度场,温度场程序写好了,但是我耦合不进去.本科毕设题目,但是专业方向不是这个,自己弄的一直在报错,求助大佬帮帮,可以私我拿程序.
- 格子玻尔兹曼方法（LBM）的学习笔记1（附Couette流源代码及解析）
笔记目录 关于学习的教材及说明 在学习之前大致将流体力学学了一下包括一些概念的理解和重要的公式,在看这本<The Lattice Boltzmann Method Principles and ...
最新文章
- 傅里叶(FFT)+小波变换+数据压缩
- Python实现信息自动配对爬虫排版程序（附下载）
- IC/FPGA笔试/面试题分析（十一）基础概念（三态门等）
- 吴裕雄 python 机器学习——数据预处理标准化StandardScaler模型
- optee中的arm64的virt_to_phys的实现
- python字符串与列表的相互转换
- php 图片上传预览(转)
- vs最好的版本_Win10 环境下，LightGBM GPU 版本的安装
- python安装不了jupyter_python学习笔记——Windowns下Python3之安装jupyter
- Xadmin添加用户小组件出错
- Kaggle入门——房价预测
- python画图xlable显示中文_xlabel和ylabel超出绘图区域，无法在figu中完全显示
- java 记账系统_案例分享用java开发实现一个记账系统（代码全）
- matlab将函数展开成幂级数,解析函数展开成幂级数的方法分析.doc
- 【PHP练习】每日词汇，随机产生10个单词，方便备考随时背诵(php+html+css)
- 计算机丢失msvcp90dll怎么办,msvcp90.dll
- Crust Network 与京湘豫等地区块链名企、投资人考察广西区块链科创园
- 弗吉尼亚大学计算机就业如何,假设你是新华中学的学生李华，高中毕业后想到美国弗吉尼亚大学(University of Virginia)计算机专业深造...
- Gos —— shell程序
- DAEFRHDSGYEVHHQKLVFFAEDV|138648-77-8