加重力的方腔流！
%还是老规矩先宣传一下QQ群： 格子玻尔兹曼救星：293267908。
%加了重力之后上侧的密度会有逐渐的减轻，因为LBM中，压力和密度是耦合的！
% ^y
% 7 3 6
% \|/
% 4-1-2 ---> x
% /|\
% 8 5 9
%
%
% INPUTS: rho The fluid density, which is a scalar.
% dx, dy, The spatial step-size.
% Lx, Ly The length of the domain.
% OUTPUTS:
%
% Authors: Qiujie Meng, Copyright 2019-2021
clear all
clc
% GENERAL FLOW CONSTANTS
Lx= 100;%Lx= 0~x_max
Ly= 100;
m=100; % m= 1~x_max
n=100;
nu= 0.01; % kinematic viscosity
dt=1;
dx=1;
dy=1;
gravity=5*10^(-6);
VSSum2=[];
x=1:dx:100;
y=1:dy:100;
rhoo=5;
u0 = 0.1 ; %0.1 % maximum velocity of Poiseuille inflow
Re = u0*Lx/nu ; % Reynolds number
omega = 1. / (3*nu+1./2.); % relaxation parameter
maxT = 40000 ; % total number of iterations
tPlot = 5; % cycles
cs=1/(3)^0.5;
% D2Q9 LATTICE CONSTANTS
wt = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cx = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
cy = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
%存一下x坐标数据 x coordinate data
data_x1 = fopen('D:\MyStudy\MATLAB\MatIB-master\tests\x1.txt', 'w');
%存一下y坐标数据 y coordinate data
data_y1 = fopen('D:\MyStudy\MATLAB\MatIB-master\tests\y1.txt', 'w');
%存一下x方向速度数据. x coordinate velocity data
data_velocity_u1 = fopen('D:\MyStudy\MATLAB\MatIB-master\tests\velocity_u1.txt', 'w');
%存一下y方向速度数据.
data_velocity_v1 = fopen('D:\MyStudy\MATLAB\MatIB-master\tests\velocity_v1.txt', 'w');
% INITIAL CONDITION: (rho=0, u=0) ==> fIn(i) = t(i)
%fIn = reshape( wt' * ones(1,Lx*Ly), 9, Lx, Ly);
rho=ones(Lx,Ly).*rhoo;
for i=1:9
fIn(i,:,:)=wt(i).*rho;
end
ux=zeros(Lx,Ly-1);
uy=zeros(Lx,Ly);
ux(:,Ly)=u0;
% MAIN LOOP (TIME CYCLES)---------------------
% MAIN LOOP (TIME CYCLES)---------------------
% MAIN LOOP (TIME CYCLES)---------------------
tic
for cycle = 1:maxT
mm=max(max(rho))%观察最大密度
if mm>100%判断仿真是否发散
break;
end
% 1、COLLISION STEP
% % 方法1
t1 = ux.^2+uy.^2;
for i = 1:9
t2 = ux.*cx(i) + uy.*cy(i);
feq(i,:,:)=wt(i).*rho.*(1+3*t2+4.5*t2.^2-1.5*t1);
force_latt(i,:,:) = (1- 0.5 * omega)*wt(i)* (3 * ((cy(i)-uy).* gravity) + 9 * (cx(i)*ux+ cy(i)*uy) * (cy(i)*gravity));
Temp(i,:,:)=(1-omega).*fIn(i,:,:)+omega.*feq(i,:,:)+force_latt(i,:,:);
end
p1=fIn; fIn =Temp; Temp = p1;
% % 方法2
% t1 = ux.^2+uy.^2;
% for i = 1:9
% t2 = ux.*cx(i) + uy.*cy(i);
% feq(i,:,:)=wt(i).*rho.*(1+3*t2+4.5*t2.^2-1.5*t1);
% Temp(i,:,:)=(1-omega).*fIn(i,:,:)+omega.*feq(i,:,:);
% end
% p1=fIn; fIn =Temp; Temp = p1;
% 2、STREAMING STEP--by Meng--莫哈默德LBM计算方法
% % Way 1
% for j1=1:Ly
% for i1= Lx:-1:2
% % Right->Left
% fIn(2, i1, j1) =fIn(2,i1-1, j1) ;
% end
% for i2= 1:(Lx-1)
% % Left->Right
% fIn(4,i2, j1) =fIn(4,i2+1, j1) ;
% end
% end
% % Top->Bottom
% for j1=Ly:-1:2
% for i1= 1:1:Lx
% fIn(3,i1,j1) =fIn(3,i1,j1-1) ;
% end
% for i2= Lx:-1:2
% fIn(6,i2,j1) =fIn(6,i2-1,j1-1) ;
% end
% for i3= 1:1:Lx-1
% fIn(7,i3,j1) =fIn(7,i3+1,j1-1) ;
% end
% end
% % Bottom->Top
% for j1=1:1:Ly-1
% for i1= 1:1:Lx
% fIn(5,i1,j1) =fIn(5,i1,j1+1) ;
% end
% for i2= 1:1:Lx-1
% fIn(8,i2,j1) =fIn(8,i2+1 ,j1+1) ;
% end
% for i3= Lx:-1:2
% fIn(9,i3,j1) =fIn(9,i3-1 ,j1+1) ;
% end
% end
% Way 2
% fIn(2,2:Lx,:)=fIn(2,1:Lx-1,:);
% fIn(3,:,2:Ly)=fIn(3,:,1:Ly-1);
% fIn(4,1:Lx-1,:)=fIn(4,2:Lx,:);
% fIn(5,:,1:Ly-1)=fIn(5,:,2:Ly);
% fIn(6,2:Lx,2:Ly)=fIn(6,1:Lx-1,1:Ly-1);
% fIn(7,1:Lx-1,2:Ly)=fIn(7,2:Lx,1:Ly-1);
% fIn(8,1:Lx-1,1:Ly-1)=fIn(8,2:Lx,2:Ly);
% fIn(9,2:Lx,1:Ly-1)=fIn(9,1:Lx-1,2:Ly);
%Way 3
for i=1:9
fIn(i,:,:) = circshift(fIn(i,:,:), [0,cx(i),cy(i)]);
end
% 3、BOUNDARY CONDITIONS
for j=1:Ly
% West boundary
fIn(2,1,j)=fIn(4,1,j);
fIn(6,1,j)=fIn(8,1,j);
fIn(9,1,j)=fIn(7,1,j);
% East boundary
fIn(4,Lx,j)=fIn(2,Lx,j);
fIn(8,Lx,j)=fIn(6,Lx,j);
fIn(7,Lx,j)=fIn(9,Lx,j);
end
% South boundary
for i=1:Lx
fIn(3,i,1)=fIn(5,i,1);
fIn(6,i,1)=fIn(8,i,1);
fIn(7,i,1)=fIn(9,i,1);
% Moving lid--north boundary
rhon=fIn(1,i,Ly)+fIn(2,i,Ly)+fIn(4,i,Ly)+2.*(fIn(3,i,Ly)+fIn(7,i,Ly)+fIn(6,i,Ly));
fIn(5,i,Ly)=fIn(3,i,Ly);
fIn(9,i,Ly)=fIn(7,i,Ly) +rhon*u0/6;
fIn(8,i,Ly)=fIn(6,i,Ly) -rhon*u0/6;
end
% 4、 MACROSCOPIC VARIABLES
for i=1:Lx
for j=1:Ly
ssum1=0;
for i0=1:9
ssum1 = ssum1+fIn(i0,i,j);
end
rho(i,j)=ssum1;
end
end
rho(:,Ly)=fIn(1,:,Ly)+fIn(2,:,Ly)+fIn(4,:,Ly)+2*(fIn(3,:,Ly)+fIn(7,:,Ly)+fIn(6,:,Ly));
ux = reshape ( (cx * reshape(fIn,9,Lx*Ly)),1, Lx,Ly); %
uy = reshape ( (cy * reshape(fIn,9,Lx*Ly)),1, Lx,Ly);
%rho = sum(f_in);
%ux = reshape ( (cx_f * reshape(f_in,9,lx*ly)), 1,lx,ly) ./rho;
%uy = reshape ( (cy_f * reshape(f_in,9,lx*ly)), 1,lx,ly) ./rho;
ux=squeeze(ux)./rho; % Reduce a matrix dimension: imporant!!
uy=squeeze(uy)./rho; % Reduce a matrix dimension: imporant!!
% 5、VISUALIZATION
%if (mod(cycle,tPlot)==0)
u = reshape(sqrt(ux.^2+uy.^2),Lx,Ly);
imagesc(u');
axis equal off; drawnow
%end
%6、 每 saveStep 步保存下误差，检查收敛性
temp1=0;
for i=1:Lx
for j=1:Ly
temp1=temp1+sqrt(ux(i,j)^2+uy(i,j)^2); %节点的速度和
end
end
VSSum2=[VSSum2;temp1];
%7、速度和、质量和
uv=sqrt(ux.^2+uy.^2);
TmpMass = 0;
Den2=sum(sum(fIn(1,:,:)+fIn(2,:,:)+fIn(3,:,:)+fIn(4,:,:)+fIn(5,:,:)+fIn(6,:,:)+fIn(7,:,:)+fIn(8,:,:)+fIn(9,:,:)));
TmpMass=TmpMass + Den2;
end
toc
% 6、 Consequence
for i=1:Lx
for j=1:Ly
fprintf(data_x1,'%f\n',i); %存x。
fprintf(data_y1,'%d\n',j); %存y。
fprintf(data_velocity_u1,'%f\n', ux(i,j));
fprintf(data_velocity_v1,'%f\n', uy(i,j));
end
end
figure
u11=uy(:,50)./u0; %y轴中间
u11=u11';
plot(x./101,u11,'*-')
hold on
u22=ux(50,:)./u0;
u22=u22';
plot(y/101,u22,'o-')
% u3=u(:,300)/uo;
% plot(y,u3,'--')
hold off
figure
uv=sqrt(ux.^2+uy.^2);
contourf(1:Ly,1:Lx,uv)
figure
quiver(ux,uy)