SIMPLE算法Matlab代码

%--------------------------------------------------------% % SIMPLE Algorithm - staggered mesh % Incompressible Lid-driven Flow in Cavity %--------------------------------------------------------%

%--------------------------------------------------------%
%       SIMPLE Algorithm - staggered mesh
%       Incompressible Lid-driven Flow in Cavity
%--------------------------------------------------------%
close all
clear
clc
disp（‘SIMPLE Algorithm‘）;

%% Setting
n=51;
m=n;
dx=1/（n-1）;
dy=1/（m-1）;
Re=100;
dt=0.001;
tmax=10000;
k1max=300;
k2max=300;
alpha=1;
alpha_p=1;
relax=1
r_p_critial=1e-2;
r_u_critial=1e-2;
r_v_critial=1e-2;
r_pc_critial=1e-2;

%% Initialization
u=zeros（n-1m）;
u（:m）=1;      % lid
v=zeros（nm-1）;
p=zeros（nm）;
r_p=zeros（tmax1）;
r_u=zeros（tmax1）;
r_v=zeros（tmax1）;

%% Loops
tic
fprintf（‘Iter\tP-residual\t\tU-residual\t\tV-residual\n‘）;
A1=1/dx^2;
A2=1/dy^2;
A3=1/（2*A1+2*A2）;
A=1/dt+2/Re*（A1+A2）;
for k=1:tmax
    %% Compute tentative velocities uv
    u1=zeros（n-1m）;
    u1（:m）=1;
    v1=zeros（nm-1）;
    for k1=1:k1max
        for i=2:n-2
            for j=2:m-1
                vn=（v（ij）+v（i+1j））/2;
                vs=（v（ij-1）+v（i+1j-1））/2;
                a1=-u（i+1j）/（2*dx）+A1/Re;
                a2=u（i-1j）/（2*dx）+A1/Re;
                a3=-vn/（2*dy）+A2/Re;
                a4=vs/（2*dy）+A2/Re;
                u1（ij）=（1-relax）*u1（ij）+relax*（u（ij）/dt+a1*u1（i+1j）+a2*u1（i-1j）+a3*u1（ij+1）+a4*u1（ij-1）-（p（i+1j）-p（ij））/dx）/A;
            end
        end
        r_u1=norm（u-u1）;
        if r_u1    end
    for k1=1:k1max
        for i=2:n-1
            for j=2:m-2
                ue=（u（ij）+u（ij+1））/2;
                uw=（u（i-1j）+u（i-1j+1））/2;
                b1=-v（ij+1）/（2*dy）+A2/Re;
                b2=v（ij-1）/（2*dy）+A2/Re;
                b3=-ue/（2*dx）+A1/Re;
                b4=uw/（2*dx）+A1/Re;
                v1（ij）=（1-relax）*v1（ij）+relax*（v（ij）/dt+b1*v1（ij+1）+b2*v1（ij-1）+b3*v1（i+1j）+b4*v1（i-1j）-（p（ij+1）-p（ij））/dy）/A;
            end
        end
        r_v1=norm（v-v1