wave_vorticity.m

clear 
clc
RE=5000;
niu=1.139*10^（-6）;
TN=500000;
err=0.0001;
rho=0.05;
L=1;
h=0.01;
N=L/h+1;
u_0=niu*RE/L;
Tw=5000;
E1=0;
E2=1;
counter=0;
U=zeros（NN）;%流量
U0=zeros（NN）;
Z=zeros（NN）;%涡量
Z0=zeros（NN）;
Vx=zeros（NN）;
Vy=zeros（NN）;

%t=0时刻初始值即U Z VX VY均为0;
for t=1:TN
    %求边界上的涡量值
    %DC
    for i=2:N-1
        Z（iN）=2*（（U（iN）-U（iN-1）-u_0*h）/h^2）;
    end
    %AB
    for i=2:N-1
        Z（i1）=2*（（U（i1）-U（i2））/h^2）;
    end
    %AD
    for j=2:N-1
        Z（1j）=2*（（U（1j）-U（2j））/h^2）;
    end
    %BC
    for j=2:N-1
        Z（Nj）=2*（（U（Nj）-U（N-1j））/h^2）;
    end
    %四个顶点的涡量值
    Z（11）=0.5*（Z（12）+Z（21））;
    Z（1N）=0.5*（Z（1N-1）+Z（2N））;
    Z（N1）=0.5*（Z（N-11）+Z（N2））;
    Z（NN）=0.5*（Z（NN-1）+Z（N-1N））;
    
    %由2.7.24求内点的新涡量值
    for i=2:N-1
        for j=2:N-1
            u=max（max（abs（Vx）））;
            v=max（max（abs（Vy）））;
            dt=min（h^2/4/niuh/（u+v））;
            Z（ij）=Z（ij）+（Vx（ij）*（Z（i-1j）-Z（i+1j））+Vy（ij）*（Z（ij-1）-Z（ij+1）））...
                *（dt/2/h）+（niu*dt/h^2）*（（Z（i+1j）-2*Z（ij）+Z（i-1j））+（Z（ij+1）-2*Z（ij）