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

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

 March 3, 2016, 22:46 Simple stream function vorticity formulation for flow past a rectangle [Matlab] #1 New Member   Kohbodi Join Date: Jun 2015 Posts: 6 Rep Power: 3 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;

 March 5, 2016, 04:44 #2 Senior Member     Join Date: May 2012 Posts: 156 Rep Power: 6 Hey, I have not tested your code, but I wonder if you have implemented the vorticity at the protruding corners correctly.

 March 5, 2016, 05:04 #3 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 2,419 Rep Power: 31 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

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post cain Main CFD Forum 5 July 19, 2012 10:05 ivanyao OpenFOAM Running, Solving & CFD 6 September 5, 2008 20:50 Christian Main CFD Forum 2 February 27, 2007 07:27 Hong Main CFD Forum 1 February 8, 2000 18:18 Willard Lee Main CFD Forum 1 January 20, 1999 12:03

All times are GMT -4. The time now is 06:09.