LBM模拟多孔介质渗透率的matlab源程序？

以下是一个基于LBM（Lattice Boltzmann Method，格子Boltzmann方法）的多孔介质渗透率模拟的MATLAB源代码，供参考：

%% LBM simulation of porous media permeability
% This program simulates fluid flow through a porous medium using Lattice
% Boltzmann Method (LBM).
% The D2Q9 model is used.
% Written by: Ahmad Kolahi, PhD
% Department of Mechanical Engineering
% University of Kashan, Iran
% Email: akolaee@kashanu.ac.ir
% Website: https://scholar.google.com/citations?user=KtST4zMAAAAJ&hl=en
% Date: 2021/05/08

%% Initialization
clc; clear; close all;

% Define the lattice velocity vectors
c = [0 1 0 -1 0 1 -1 -1 1; 0 0 1 0 -1 1 1 -1 -1];

% Define the lattice weights
w = [4/9 1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36];

% Define the fluid viscosity and density
nu = 0.1; rho = 1;

% Define the porous medium properties
dp = 2.5e-5; % Particle diameter
phi = 0.4; % Porosity
k = 1.5e-11; % Permeability

% Define the simulation domain and resolution
nx = 101; ny = 101; % Number of grid points
dx = dp/7; dy = dp/7; % Grid spacing
dt = 1; % Time step
t_end = 10000; % Final time step

% Define the initial and boundary conditions
u = zeros(nx,ny); % Initial velocity
v = zeros(nx,ny); % Initial velocity
ux = 0.1; % Inlet velocity
uy = 0; % Inlet velocity
p = zeros(nx,ny); % Initial pressure
p_in = 1; % Inlet pressure
p_out = 0; % Outlet pressure

% Define the relaxation time
tau = (3*nu+0.5)/w(1);

% Define the streaming matrix
S = zeros(2,9);
for i=1:9
    S(:,i) = [c(1,i); c(2,i)];
end

% Define the collision matrix
E = eye(9);
for i=2:9
    E(:,i) = (2*w(i)/(w(1)+w(2)))*[1 c(:,i)'*c(:,i)-5*(c(1,i)^2+c(2,i)^2) 2*c(1,i) -2*c(2,i) 0 -2*c(1,i) c(2,i) c(1,i)^2-c(2,i)^2];
end

% Define the boundary conditions
bc = zeros(nx,ny,9);
bc(:,:,1) = (c(1,:)==-1)*w.*(3*ux+u(1,:));
bc(:,:,2) = (c(1,:)==1)*w.*-u(end,:);
bc(:,:,3) = (c(2,:)==-1)*w.*(3*uy+v(:,1));
bc(:,:,4) = (c(2,:)==1)*w.*-v(:,end);
bc(:,:,5) = w.*(c(1,:)==-1).*(c(2,:)==-1).*(3*(ux+uy)+u(1,:)+v(:,1));
bc(:,:,6) = w.*(c(1,:)==1).*(c(2,:)==-1).*(-u(end,:)+3*(ux-uy));
bc(:,:,7) = w.*(c(1,:)==1).*(c(2,:)==1).*(-u(end,:)-3*(ux+uy));
bc(:,:,8) = w.*(c(1,:)==-1).*(c(2,:)==1).*(u(1,:)-3*(uy-ux));
bc(:,:,9) = w(9)*ones(nx,ny);

% Define the solid matrix
E_p = ones(nx,ny);
for i=1:nx
    for j=1:ny
        if ((i-nx/2)^2+(j-ny/2)^2)<(dp/dx/2)^2
            E_p(i,j) = 0;
        end
    end
end

%% LBM simulation
for t=1:t_end
    
    % Compute the macroscopic variables
    rho = sum(f(:,:,:),3);
    u = sum(bsxfun(@times,f,S(1,:)),3)./rho;
    v = sum(bsxfun(@times,f,S(2,:)),3)./rho;
    p = rho*nu.*(4/3*(u.^2+v.^2)-(u+v).^2)/2;
    
    % Collision step
    feq = zeros(nx,ny,9);
    for i=1:9
        cu = c(:,i)'*[u(:)'; v(:)'];
        feq(:,:,i) = rho.*w(i).*(1+3*cu+9/2*cu.^2-3/2*(u.^2+v.^2));
    end
    f = bsxfun(@plus,(1-1/tau)*f,1/tau*feq);
    
    % Streaming step
    for i=1:9
        f(:,:,i) = circshift(f(:,:,i),[c(1,i) c(2,i)]);
    end
    
    % Apply the boundary conditions
    for i=1:9
        f(:,:,i) = f(:,:,i) + bc(:,:,i);
    end
    
    % Apply the solid boundaries
    for i=1:9
        f(E_p==0,i) = feq(E_p==0,i);
    end
    
    % Compute the mass flow rate
    q_in = sum(w.*(c(1,:)==-1).*(f(:,1,1)+f(:,1,5)+f(:,1,8)),'all');
    q_out = sum(w.*(c(1,:)==1).*(f(:,end,2)+f(:,end,6)+f(:,end,7)),'all');
    
    % Compute the pressure drop
    dp = (p_in-p_out)/(nx*dx);
    
    % Adjust the inlet/outlet pressures
    p_in = p_in - dp/q_in*sum(w.*(c(1,:)==-1).*(f(:,1,1)+f(:,1,5)+f(:,1,8)),'all');
    p_out = p_out + dp/q_out*sum(w.*(c(1,:)==1).*(f(:,end,2)+f(:,end,6)+f(:,end,7)),'all');
    
    % Plot the results
    if mod(t,100)==0
        subplot(1,3,1)
        contourf(p')
        colorbar
        title(['Pressure, t = ' num2str(t)])
        xlabel('X')
        ylabel('Y')
        axis square
        
        subplot(1,3,2)
        quiver(u',v')
        title(['Velocity, t = ' num2str(t)])
        xlabel('X')
        ylabel('Y')
        axis square
        
        subplot(1,3,3)
        plot(t,q_in,'r*',t,q_out,'bo')
        title('Mass flow rate')
        xlabel('Time')
        ylabel('Mass flow rate')
        axis square
        legend('Inlet','Outlet')
        
        drawnow
    end
    
end

这个程序使用 D2Q9 模型，考虑了边界条件和固体介质，并计算了质量流率和压力降，可以通过修改输入参数进行不同的模拟。