CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Simple stream function vorticity formulation for flow past a rectangle [Matlab]

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   March 3, 2016, 22:46
Default Simple stream function vorticity formulation for flow past a rectangle [Matlab]
  #1
a99
New Member
 
Kohbodi
Join Date: Jun 2015
Posts: 6
Rep Power: 3
a99 is on a distinguished road
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;
a99 is offline   Reply With Quote

Old   March 5, 2016, 04:44
Default
  #2
Senior Member
 
Simbelmynė's Avatar
 
Join Date: May 2012
Posts: 157
Rep Power: 6
Simbelmynė is on a distinguished road
Hey,

I have not tested your code, but I wonder if you have implemented the vorticity at the protruding corners correctly.
Simbelmynė is offline   Reply With Quote

Old   March 5, 2016, 05:04
Default
  #3
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,700
Rep Power: 33
FMDenaro will become famous soon enoughFMDenaro will become famous soon enough
I think no one will check you code line-by-line, therefore if you problem is a numerical instability check if your time step is in the stability region.

Furthermore, I suggest to perform one single time step a plotting vorticity and stream function
FMDenaro is online now   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Streamfunction? cain Main CFD Forum 5 July 19, 2012 10:05
Problem with compile the setParabolicInlet ivanyao OpenFOAM Running, Solving & CFD 6 September 5, 2008 20:50
Droplet Evaporation Christian Main CFD Forum 2 February 27, 2007 07:27
Stream Function Formulation Hong Main CFD Forum 1 February 8, 2000 18:18
Stream function and vorticity for turbulence flow? Willard Lee Main CFD Forum 1 January 20, 1999 12:03


All times are GMT -4. The time now is 16:39.