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

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.