[JungleCrow]

Hey,I am a beginner in Matlab and recently I wrote a code for finding flow around a square Cross-section using Vorticity Stream Function.This is how I went about:
1:Assumed value of Stream Function on square Boundary as 0.
2:Found out the Stream Function value on top and Bottom Boundary of domain by finding the volume flow using Inlet Velocity and length of inlet Boundary.
3:Assumed linear variation of Stream Function on inlet Boundary as velcity is constant.
4:Vorticity is assumed to be 0 on left(inlet),top & bottom Boundary[assumed domail large enough such that effect of cylinder not there)
5:Found out Boundary conditions for Vorticity on 4 sides of square.(similar to Lid Driven Cavity) problem.
6:Used gauss Seidal method with over-relaxation of 1.5 for faster Convergence.

But unfortunately i am unable to get proper results even after running the code for hours. I will be highly grateful if someone could tell me where I am going wrong and suggest right way to proceed.
Thank you in advance.

grid points for square is (301,301) to (331,331)
Re=500


% Date - 09/07/2015

clear all
clc

%% Input Data
dimr=0.3; % Dimension of the riser
cv=0.0017; % Current velocity
kv=1.004e-6; % Kinematic viscosity
Re=500; % Reynolds number

nx=1831; % Grid dimension in X-direction
ny=631; % Grid dimension in Y-direction
x=0:0.01:18.3; % Grid size in X-direction
y=0:0.01:6.3; % Grid size in Y-direction
ni=100; % Number of iterations
nt=100;
err=0.01; % Acceptable Error
%Forcesw=0;
%Forceaw=0;
x2d=zeros(ny,nx);
y2d=zeros(ny,nx);
t=1:0.01:1000;
dt=t(2)-t(1);
ds=0.1;
SOR=1.5;%Successive Over Relaxation
vort=zeros(ny,nx);
U1=.0017;%Inlet Velocity
pn=zeros(ny,nx);
phi(ny,1:nx)=3;% Stream function Value at Top Boundary
phi(1,1:nx)=-3;% Stream Function Value at Bottom Boundary
for a=1:ny% Assuming Linear Variation in Stream function value for left Boundary
phi(a,1)=0.952*y(a)-3;
end
for z=301:331
for b=301:331
phi(z,b)=0;
end
end




%% Calculations
% Stream Function calculation
for ts=1:nt % time steps start
for iter=1:ni % stream func. calculation
K=phi; % Successive Over-Relaxation
for i=2:ny-1;
for j=2:nx-1;
phi(i,j)=0.25*SOR*(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1)+ds*ds*vort(i,j))+(1-SOR)*phi(i,j);
end
end
error=0.0;
for i=1:ny
for j=1:nx
error=error+abs(K(i,j)-phi(i,j));
end
end
if error<=err,
break;
end%if iteration converge terminate the loop
end

%Vorticity% making vorticity 0 and left,top & bottom boundary of
%domain
for i=1
for j=1:nx
vort(i,j)=0.0;
end
end
for i=ny
for j=1:nx
vort(i,j)=0.0;
end
end
for j=1
for i=1:ny
vort(i,j)=0.0;
end
end
for i=301:331 % Boundary condition for Vorticity on left side of square cylinder(cross-section)
for j=301
vort(i,j)=2*phi(i-1,j)/(ds*ds);
end
end
for i=301:331 % Boundary condition for Vorticity on Right side of square cylinder
for j=331
vort(i,j)=-2*phi(i+1,j)/(ds*ds);
end
end
for i=301 % Boundary condition for Vorticity on Bottom side of square cylinder
for j=301:331
vort(i,j)=2*phi(i,j-1)/(ds*ds);
end
end
for i=331 % Boundary condition for Vorticity on Bottom side of square cylinder
for j=301:331
vort(i,j)=-2*phi(i,j+1)/(ds*ds);
end
end
for j=1:nx%assuming vorticity at right boundary(column)= vorticity at 1 column before
vort(nx)=vort(nx-1);
end


%................................................. .................
for i=2:ny-1;
for j=2:nx-1 % compute
K(i,j)=-0.25*((phi(i,j+1)-phi(i,j-1))*(vort(i+1,j)-vort(i-1,j))-(phi(i+1,j)-phi(i-1,j))*(vort(i,j+1)-vort(i,j-1)))/(ds*ds)+(1/Re)*(vort(i+1,j)+vort(i-1,j)+vort(i,j+1)+vort(i,j-1)-4.0*vort(i,j))/(ds*ds);
end
end

% Updating the vorticity
vort(2:ny,2:nx)=vort(2:ny,2:nx)+dt*K(2:ny,2:nx);

%for i=1:nx
%for j=1:ny
%x2d(i,j)=x(i);
%y2d(i,j)=y(j);
%end
%end

%calculation of U and V
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>

for i=2:ny-1
for j=2:nx-1
u(i,j)=(phi(i,j+1)-phi(i,j))/(2*ds);
v(i,j)=(phi(i+1,j)-phi(i,j))/(2*ds);
u(ny,: )=U1;
v(ny,: )=0;
u(:,1)=U1;
v(:,1)=0;
u(1,: )=U1;
v(1,: )=0;
for i=301:331
for j=301:331
u(i,j)=0;
v(i,j)=0;
end
end
end
end
end

%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>>>
% Pressure calculation
h=zeros(ny,nx);
for i=2:ny-1
for j=2:nx-1
h(i,j)=(((phi(i-1,j)-2*phi(i,j)+phi(i+1,j))/(ds*ds))*((phi(i,j-1)-2*phi(i,j)+phi(i,j+1))/(ds*ds)))-(phi(i+1,j+1)-phi(i+1,j-1)-phi(i-1,j-1))/(4*(ds*ds));
p(i,j)=(0.25*(pn(i+1,j)+pn(i-1,j)+pn(i,j+1)+pn(i,j-1))-0.5*((h(i,j)*ds*ds*ds*ds)));
end
pn(i,j)=p(i,j);
end

[FMDenaro]

it will be much more useful if you post the image of the stream function after one time-step....

However, if you are simulating a bluf-body inside a 2D channel, you can fix psi=0 on the bottom wall of channel and normalize the flow rate to have psi =1 on the top channel wall. Inflow velocity for a plug-flow is a linear psi function between 0 and 1 with zero vorticity. Be carefull at the outflow conditions.

Then, the psi on the bluff-body is a value dependent on the flow.
If the geometry is exactly symmetric and the flow is steady, you can fix psi=0.5 and compute the vorticity with some formula, for example the Thom formula.
But if the flow is unsteady or the bluff body is not at the center, you have to compute dynamically the value of psi.

[JungleCrow]

Thankyou for replying.
In my case geometry is symmetric.I had taken the stream function value 0 at square just because it will be convenient to write the boundary conditions on vorticity (the the 1st term and the velocity term becomes 0 on expansion of stream function using Taylor series)

After around 5 minutes of running the error was around 2500.The program reached 3rd time step after around 45 hours of running.With error of 242. How should I modify the program so that it converges faster?Should I increase the SOR parameter? Shall I adopt a new mathematical method to solve the equation like Fourier Transformation(i used Gauss Seidal for this).

This is how the stream function looks like after 3 time steps.Error is 242.
Stream Function with X


Stream Function with Y
http://imageshack.com/a/img912/8483/eO5Hh6.jpg

Stream function with X,Y
http://imageshack.com/a/img673/1694/poYF7f.jpg

[FMDenaro]

1) in viscous flow vorticity is not zero at a wall!
2) I do not understand the geometry defined by the streamlines. Where the bluff body is?
3) outflow conditions are wrong


I suppose you have to revise many issues in your code... some theoretical errors as well as some simple bug coding can be present.

[Alex C.]

It is not clear to me if you are trying to solve 2D channel or 2D open flow around bluff body.

If you are doing channel, the correct boundary conditions is not zero vorticity at wall. The correct boundary condition would be an equation between vorticity flux at wall and pressure gradient. Refer to a fluid mechanic textbook on this subject. I have a preference for Kundu et al. textbook.

If you are doing open flow, and your domain is large enough (width to be verified when solution seems right), then zero vorticity might be fine (? Can anyone confirm) Again a textbook on fluid mechanics can help.

As previously said by prof Denaro, the outflow is wrong.

The inflow seems ok, but looks odd because the boundary condition on the bluff body are not set correctly. Reading your code, I see you set the no-slip BC in velocity variables, but you don't recalculate the vorticity and stream function from the modified velocity field after the BC is applied. If you proceed this way, the vorticity and stream functions don't "see" the square cross-section from one iteration to the other. I suspect that what we see in the stream lines you showed is related to initial condition.

When you are doing stream relaxation, the maximum number of iteration is set to 100. You might want to check how your solver exits the relaxation loop. You might end up running 100 iterations and the error is still very large. Also the error threshold might be too big (or too small).

Last comment : "phi(a,1)=0.952*y(a)-3;"
I don't like this line, as it is linked to parameters that you might want to change, but is hard-coded. What if you change the top and bottom stream function value? What if you change your grid spacing? I haven't look for other such things, but there might be others, and you should be careful about them.

The long post might be depressing and hard on the work done. I wish to encourage you to proceed. CFD is not something easy, and coding problems like this one is definitely a challenge that will help you learn a lot. Keep up the good work.

Cheers.

[FMDenaro]

I suggest tryng first a simple case: inviscid flow around the bluff-body, this way you have simply to solve the problem


Lap Psi = 0

v.n = 0 on wall
v.n = vn on inflow/outflow

a Jacobi or Gauss-Seidel algorithm is fine to solve the linear system.
When you did that, we can check better the streamline.
Again, please highlight the body position

[JungleCrow]

Thankyou for taking the pain of reading the code.I will modify the code,keeping in mind your suggestions.