1
$\begingroup$

I want to calculate the stream function $\psi$ starting from a velocity field $(u,v)$ (such that $u=-\frac{\partial\psi}{\partial y}$ and $v=\frac{\partial\psi}{\partial x}$). I thus calculate the vorticity $\omega$ (such that $\omega=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$) and then solve the Poisson equation $\nabla^2 \psi=\omega$.

To easily specify the boundary conditions, I choose a successive over-relaxation (SOR) approach.

The code seemingly runs smoothly. However, when I check if the solution is correct by differentiating $\psi$, its derivatives do not match the initial velocities $(u,v)$

Is this due to some intrinsic limitation of the SOR?

u

v

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')
$\endgroup$
  • $\begingroup$ You're asking a debugging question. If you wonder whether it's due to the SOR method, have you tried to replace that part of your code by, say, a direct solver for the system of equations? $\endgroup$ – Wolfgang Bangerth Apr 4 '18 at 2:56

This site is temporarily in read only mode and not accepting new answers.

Browse other questions tagged .