加重力的方腔流！

%还是老规矩先宣传一下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