我想从速度场（这样和）。因此，我计算涡量（使得），然后求解泊松方程。$\psi$$(u,v)$$u=-\frac{\partial\psi}{\partial y}$$v=\frac{\partial\psi}{\partial x}$$\omega$$\omega=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$$\nabla^2 \psi=\omega$

为了轻松指定边界条件，我选择了连续过松弛 (SOR) 方法。

代码看起来运行流畅。检查解决方案是否正确时不匹配$\psi$$(u,v)$

这是由于 SOR 的一些内在限制吗？

close all
clearvars
clc


%coordinates
x=1:10:600; x = repmat(x',1,50); 
y=1:10:500; y = repmat(y,60,1);
dx = 10*ones(60,50); %x spacings
dy = 5*ones(60,50); %y spacings

%define velocities
u = y.^3;
v = x.^3;
figure;pcolor(x,y, u);title('u')
figure;pcolor(x,y, v);title('v')

%vorticity
[dudy,~] = gradient(u, dy(1,1),dx(1,1));
[~,dvdx] = gradient(v, dy(1,1),dx(1,1));
figure;pcolor(x,y, dudy);title('dudy')
figure;pcolor(x,y, dvdx);title('dvdx')

omega = dvdx - dudy;

figure;pcolor(x,y, omega);title('omega');colormap(jet(10));colorbar


%iterative solver (SOR)
psi=rand(62,50); %initial guess

beta=dx./dy; %spacing ratio

%boundary conditions in y
omega(:,1:2)=0;
omega(:,end-1:end)=0;

%periodic boundary conditions in x -- append/prepend values on the other side
omega = [omega(end-1:end,:); omega; omega(1:2,:)]; 
beta = [beta(end-1:end,:); beta; beta(1:2,:)];

[M,N]=size(psi);

w=1.4; %relaxation factor

figure
for iter=1:1000
  for ii=3:M-2
    for jj=3:N-2

%%fourth order formula 
psi(ii,jj) = (1-w)*psi(ii,jj) + w*...
                    ( 16*( psi(ii+1,jj) + psi(ii-1,jj) ) - psi(ii+2,jj) - psi(ii-2,jj) +... 
                    beta(ii,jj)^2*( -psi(ii,jj+2) + 16*(psi(ii,jj+1) + psi(ii,jj-1)) -psi(ii,jj-2) ) +...
                    12*dx(ii,jj)^2*omega(ii,jj) )...
                   /(30*(1+(beta(ii,jj)^2))); 

%%second order formula               
%          psi(ii,jj) = (1-w)*psi(ii,jj) + ...
%              w*( (psi(ii+1,jj)+beta(ii,jj)^2*psi(ii,jj+1)+psi(ii-1,jj)+beta(ii,jj)^2*psi(ii,jj-1)) -...
%              2*dx(ii,jj)^2*omega(ii,jj) )...
%              /(2*(1+(beta(ii,jj)^2)));

    end
  end

      if mod(iter,10)==0
        pcolor(x,y, psi(2:end-1,:))
        colorbar
        title(num2str(iter))
        drawnow
      end

end
psi = psi(2:end-1,:); %resize back

figure;pcolor(x,y, psi);title('psi')


%%%test
[dpsidy,dpsidx] = gradient(u, dy(1,1),dx(1,1));


figure;hold on
subplot(1,2,1)
pcolor(x,y, u); colorbar; title('u orig')
subplot(1,2,2)
pcolor(x,y, -dpsidy); colorbar; title('u reconstructed')

figure;hold on
subplot(1,2,1)
pcolor(x,y, v); colorbar; title('v orig')
subplot(1,2,2)
pcolor(x,y, dpsidx); colorbar; title('v reconstructed')