CFD Online Discussion Forums - Simple stream function vorticity formulation for flow past a rectangle [Matlab]

Hi,

I am working on a 2D stream function vorticity formulation to show the the streamlines and vorticity around a rectangle (a very simplified Ahmed body). Basically it is rectangle situated on the bottom of the domain.

I started working based on code I found for the lid driven cavity problem and edited to my needs.

The size of the domain was created so that the distance from the front edge of the bluff body with length L is 3L from the left wall and 7L from right wall measured from end of the body. The height of the domain is 5h while h is the height of the bluff body.

I have been editing this code for a while now and still keeps blowing up.

-Is the the BCs?
-Is it my parameters?

Any help or guidance is much appreciated. Here is the code (MaxStep is set to 11 since it blows up after this step):

clc
clear all
close all
clf;
nx=112;
ny=51;
MaxStep=11; % Number of iterations
Visc= 0.1; % Viscosity
dt= 0.001; % Time Step

% Governing parameters parameters for SOR iteration
MaxIt=100;
Beta=1.5;
MaxErr=0.001;

sf=zeros(nx,ny); % Stream function old
vt=zeros(nx,ny); % Vorticity
w=zeros(nx,ny); % Stream function new
h=1.0/(nx-1); % Step size
% dx=1.0/(nx-1); % Step size X direction
% dy=1.0/(ny-1); % Step size Y direction
t=0.0; % Initial Time

% Stream function initilization
for i = 1:nx; for j = 1:ny;
sf(i,j+1) = ((1/ny)*j);
end
end

for istep=1:MaxStep, % Start of the time integration
for iter=1:MaxIt , % Solve for the stream function by SOR iteration
sf(31:41,1:11) = 0; % Set stream function to 0 for Obstacle location
w=sf; % Update stream fucntion
for i=2:nx-1;for j=2:ny-1
sf(i,j)=0.25*Beta*(sf(i+1,j)+sf(i-1,j)...
+sf(i,j+1)+sf(i,j-1)+h*h*vt(i,j))+(1.0-Beta)*sf(i,j);
end;end;
Err=0.0;for i=1:nx;for j=1:ny, Err=Err+abs(w(i,j)-sf(i,j)) ;end;end;
if Err<= MaxErr,break,end % stop if iteration has converged
end;

sf(31:41,1:11) = 0; % Set stream function to 0 for Obstacle location
vt(31:41,1:11) = 0; % Set vorticity to 0 for Obstacle location

vt(2:nx-1,1)=-2.0*sf(2:nx-1,2)/(h*h); % vorticity on bottom wall
% vt(2:nx-1,ny)=-2.0*sf(2:nx-1,ny-1)/(h*h); % vorticity on top wall
% vt(1,2:ny-1)=-2.0*sf(2,2:ny-1)/(h*h); % vorticity on left wall
vt(nx,2:ny-1)=-2.0*sf(nx-1,2:ny-1)/(h*h); % vorticity on right wall

vt(31,1:11)=-2.0*sf(30,1:11)/(h*h); % vorticity on left block
vt(41,1:11)=-2.0*sf(42,1:11)/(h*h); % vorticity on right block
vt(31:41,11)=-2.0*sf(31:41,12)/(h*h); % vorticity on top block

%Compute the right hand side of the vorticity equation
for i=2:30;for j=2:ny-1
w(i,j)=-0.25*((sf(i,j+1)-sf(i,j-1))*(vt(i+1,j)-vt(i-1,j))...
-(sf(i+1,j)-sf(i-1,j ))*(vt(i,j+1)-vt(i,j-1)))/(h*h)...
+Visc*(vt(i+1,j)+vt(i-1,j)+vt(i,j+1)+vt(i,j-1)-4.0*vt(i,j))/(h*h);
end;end;

for i=31:41;for j=12:ny-1
w(i,j)=-0.25*((sf(i,j+1)-sf(i,j-1))*(vt(i+1,j)-vt(i-1,j))...
-(sf(i+1,j)-sf(i-1,j ))*(vt(i,j+1)-vt(i,j-1)))/(h*h)...
+Visc*(vt(i+1,j)+vt(i-1,j)+vt(i,j+1)+vt(i,j-1)-4.0*vt(i,j))/(h*h);
end;end;

for i=42:nx-1;for j=2:ny-1
w(i,j)=-0.25*((sf(i,j+1)-sf(i,j-1))*(vt(i+1,j)-vt(i-1,j))...
-(sf(i+1,j)-sf(i-1,j ))*(vt(i,j+1)-vt(i,j-1)))/(h*h)...
+Visc*(vt(i+1,j)+vt(i-1,j)+vt(i,j+1)+vt(i,j-1)-4.0*vt(i,j))/(h*h);
end;end;

% update the vorticity
vt(2:nx-1,2:ny-1)=vt(2:nx-1,2:ny-1)+dt*w(2:nx-1,2:ny-1);
t=t+dt % print out t

% Plot vorticity
figure(1);
plot(121),contourf(rot90(fliplr(vt))),axis([0,nx,0,ny]);
rectangle('Position',[31 1 10 10]);
colorbar;
grid on;
axis('equal');
title('Vorticity');

% Plot stream function
figure(2);
plot(122),contourf(rot90(fliplr(sf))),axis([0,nx,0,ny]);
rectangle('Position',[31 1 10 10])
colorbar;
grid on;
axis('equal');
title('Streamlines');

pause(0.01)
end;