给我找一个基于数字岩心模拟流体流动的matlab程序？

以下是一个基于数字岩心模拟流体流动的 MATLAB 程序：

% 设置模拟参数
nx = 100; % x方向网格数
ny = 100; % y方向网格数
nt = 100; % 时间步数
dt = 0.01; % 时间步长
dx = 1; % x方向网格间距
dy = 1; % y方向网格间距
visc = 0.1; % 粘性系数
rho = 1; % 流体密度

% 初始化速度和压力场
u = zeros(nx,ny);
v = zeros(nx,ny);
p = zeros(nx,ny);

% 设置边界条件
u(:,1) = 1; % 左边界
u(:,end) = 1; % 右边界
v(1,:) = 0; % 下边界
v(end,:) = 0; % 上边界

% 开始时间循环
for t=1:nt
    % 计算速度场
    un = u;
    vn = v;
    u(2:end-1,2:end-1) = un(2:end-1,2:end-1)...
        - un(2:end-1,2:end-1).*dt./dx.*(un(2:end-1,2:end-1)-un(1:end-2,2:end-1))...
        - vn(2:end-1,2:end-1).*dt./dy.*(un(2:end-1,2:end-1)-un(2:end-1,1:end-2))...
        - dt./rho./2./dx.*(p(3:end,2:end-1)-p(1:end-2,2:end-1))...
        + visc*dt/dx^2.*(un(3:end,2:end-1)-2*un(2:end-1,2:end-1)+un(1:end-2,2:end-1))...
        + visc*dt/dy^2.*(un(2:end-1,3:end)-2*un(2:end-1,2:end-1)+un(2:end-1,1:end-2));
    v(2:end-1,2:end-1) = vn(2:end-1,2:end-1)...
        - un(2:end-1,2:end-1).*dt./dx.*(vn(2:end-1,2:end-1)-vn(1:end-2,2:end-1))...
        - vn(2:end-1,2:end-1).*dt./dy.*(vn(2:end-1,2:end-1)-vn(2:end-1,1:end-2))...
        - dt./rho./2./dy.*(p(2:end-1,3:end)-p(2:end-1,1:end-2))...
        + visc*dt/dx^2.*(vn(3:end,2:end-1)-2*vn(2:end-1,2:end-1)+vn(1:end-2,2:end-1))...
        + visc*dt/dy^2.*(vn(2:end-1,3:end)-2*vn(2:end-1,2:end-1)+vn(2:end-1,1:end-2));
    
    % 计算压力场
    pn = p;
    div = zeros(nx,ny);
    div(2:end-1,2:end-1) = -rho.*((u(2:end-1,3:end)-u(2:end-1,2:end-1))./dx...
        +(v(3:end,2:end-1)-v(2:end-1,2:end-1))./dy);
    for q=1:20
        pn(2:end-1,2:end-1) = (pn(2:end-1,3:end)+pn(2:end-1,1:end-2)...
            +pn(3:end,2:end-1)+pn(1:end-2,2:end-1)-div(2:end-1,2:end-1).*dx^2./2)./4;
        pn(:,end) = pn(:,end-1); % 右边界
        pn(:,1) = pn(:,2); % 左边界
        pn(1,:) = pn(2,:); % 下边界
        pn(end,:) = 0; % 上边界
    end
    p = pn;
end

% 绘制速度场
[X,Y] = meshgrid(1:ny,1:nx);
figure;
quiver(X,Y,u',v');

这个程序模拟了在一个 $100\times100$ 的网格上流体的流动，使用了显式有限差分方法来求解 Navier-Stokes 方程。这个程序使用了压力-速度耦合方法来保持速度场的无散性。程序中使用了边界条件，左右边界为固定速度，上下边界为自由边界。最后，程序绘制了速度场的箭头图。