2-fold Question. Why does my code blow up/how to fix, can this be made faster? 2-D Fluid Flow around cylinder

3 views (last 30 days)
I've been staring at this code for the past week and have no idea how to fix this. I am simulating 2-D flow past a cylinder. So in my viscous solution, I need Karman vortex sheets and so forth. My omega (vorticity) calculation and other stuff ends up making my code unstable and everything becomes NaN.
I realized that there is some issue with the vorticity at the wall of the cylinder with some values becoming negative. But I don't know why or how to fix this.
I'm more interested in making the code work than efficiency, but If you have ideas how to make the code faster by vectorizing or profiling it would be a big help too. My movie code doesnt work as well either. This is a big question, I know. But help is appreciated.
clear;clc;close all
tic
cylgrid = ones(201,502);
init = 100:-1:-100;
i=1:201;
for j=1:502
cylgrid(:,j)=init;
end
U =1; %m/s
nu = 0.1544; %cm^2/sec
h = 0.01;%grid spacing .01cm
Reh = U*h/nu;%must be less than 10
% initialize matrix to keep track of indices that should be reset to zero
indcsToReset = cast(zeros(size(cylgrid)),'logical');
%initialize circle
j=500;
for ii = -100:100 %fluctuate i
for jj= 1:502 %fluctuate j
if (ii.^2 + (jj - j ./ 5).^2) < 625
cylgrid(ii+101,jj) = 0;%psi = 0 at surface and inside cylinder
% reset this index
indcsToReset(ii+101,jj) = true;
end
end
end
psi = cylgrid;%just for notation
psinew= cylgrid;
% get the number of rows and columns in psi (r==201,c==502)
[r,c] = size(psi);
error = 0.001;
tpot=0;dt=1;
while error < max(max(psinew-psi)) || tpot==0
psi=psinew;
psinew(2:r-1,2:c-2) = (1/4)*(psinew(1:r-2,2:c-2)+psi(3:r,2:c-2)+psi(2:r-1,3:c-1)+psinew(2:r-1,1:c-3));
%making sure we begin with our circle again...
psinew(indcsToReset)=0;
tpot=tpot+dt;
end
%PUT THIS BACK IN LATER
% figure(1); title('Steamline Plot');
% contourf(psinew,100);set(gca,'DataAspectRatio',[1 1 1]);
%Now for viscous flow solution
%compute the vorticity at all solid walls (the cylinder) and in the bulk
%initialize bulk vorticity
omega = zeros(201,502);
%initialize vorticity at the wall
for ii = -100:100 %fluctuate i
for jj= 1:502 %fluctuate j
if (ii.^2 + (jj - j ./ 5).^2) < 625
omega(ii+101,jj)= -(2/(h^2))*(psinew(ii+1+101,jj)+psinew(ii-1+101,jj)+psinew(ii+101,jj+1)+psinew(ii+101,jj-1)-4*psinew(ii+101,jj));
end
end
end
%finding u and v
u = psinew; v=psinew; %so that BC are in place
v(1:r-1,1:c-1) = (1/(2*h)).*(psinew(1:r-1,2:c)- psinew(1:r-1,1:c-1)); v(:,c)=v(:,c-1); v(r,:)= v(r-1,:);
u(1:r-1,1:c-1) = (1/(2*h)).*(psinew(1:r-1,1:c-1) - psinew(2:r,1:c-1));u(:,c)=u(:,c-1);u(r,:)=u(r-1,:);
%finite differencing of vorticity eqn: omega(n+1) = omega(n) +
%dt(-lapl(u*omega(n))/h - lapl(v*omega(n))/h +nu*lapl(omega(n)))
laplomega=omega;
laplomega(2:r-1,2:c-2)= (1/4)*(omega(1:r-2,2:c-2)+omega(3:r,2:c-2)+omega(2:r-1,3:c-1)+omega(2:r-1,1:c-3)-4*omega(2:r-1,2:c-2));
lapluom=zeros(r,c);laplvom=lapluom;
%upwind differencing
for ii = 2:r-1
for jj = 2:c-1
if u(ii,jj)<0
lapluom(ii,jj) = u(ii+1,j)*omega(ii+1,j)-u(ii,jj)*omega(ii,jj);
else
lapluom(ii,jj) = u(ii,jj)*omega(ii,jj)-u(ii-1,jj)*omega(ii-1,jj);
end
if v(ii,jj)<0
laplvom(ii,jj) = v(ii,jj+1)*omega(ii,jj+1) - v(ii,jj)*omega(ii,jj);
else
laplvom(ii,jj) = v(ii,jj)*omega(ii,jj) - v(ii,jj-1)*omega(ii,jj-1);
end
end
end
% lapluom(:,c) = lapluom(:,1);laplvom(:,c)=laplvom(:,1);
% lapluom(r,:)= lapluom(1,:); laplvom(r,:)=laplvom(r,:);
%omega BC, outflow-same as clmn before, top/bot-omega=0
%psi BC, outflow psinew(i,j) = 2*psi(i-1,j)-psi(i-2,j), top/bot psi = const
treal=0;dt=.0001; timestep=0;
omeganew = omega+dt*( -lapluom/h -laplvom/h +nu*laplomega);
F=1;
while error < max(max((F/4)*(psinew(1:r-2,2:c-2)+psi(3:r,2:c-2)+psi(2:r-1,3:c-1)+psinew(2:r-1,1:c-3)+(h^2).*omeganew(2:r-1,2:c-2)-4*psi(2:r-1,2:c-2)))) || treal==0
psinew(2:r-1,2:c-2) = psi(2:r-1,2:c-2)+(F/4)*(psinew(1:r-2,2:c-2)+psi(3:r,2:c-2)+psi(2:r-1,3:c-1)+psinew(2:r-1,1:c-3)+(h^2).*omeganew(2:r-1,2:c-2)-4*psi(2:r-1,2:c-2));
psinew(:,c) = 2*psi(:,c-1)-psi(:,c-2);%--OUTFLOW BC
psinew(indcsToReset)=0;
psi=psinew;
%initialize vorticity at the wall
for ii = -100:100 %fluctuate i
for jj= 1:502 %fluctuate j
if (ii.^2 + (jj - j ./ 5).^2) < 625
omega(ii+101,jj)= -(2/(h^2))*(psinew(ii+1+101,jj)+psinew(ii-1+101,jj)+psinew(ii+101,jj+1)+psinew(ii+101,jj-1)-4*psinew(ii+101,jj));
end
end
end
%now find omeganew again
u(1:r-1,1:c-1) = (1/(2*h))*(psinew(1:r-1,2:c)- psinew(1:r-1,1:c-1)); u(:,c)=u(:,c-1); u(r,:)= u(r-1,:);
v(1:r-1,1:c-1) = (1/(2*h))*(psinew(1:r-1,1:c-1) - psinew(2:r,1:c-1));v(:,c)=v(:,c-1);v(r,:)=v(r-1,:);
laplomega=omeganew; omega = omeganew;
laplomega(2:r-1,2:c-2)= (1/(h^2))*(omega(1:r-2,2:c-2)+omega(3:r,2:c-2)+omega(2:r-1,3:c-1)+omega(2:r-1,1:c-3)-4*omega(2:r-1,2:c-2));
for ii = 2:r-1
for jj = 2:c-1
if u(ii,jj)<0
lapluom(ii,jj) = u(ii+1,j)*omega(ii+1,j)-u(ii,jj)*omega(ii,jj);
else
lapluom(ii,jj) = u(ii,jj)*omega(ii,jj)-u(ii-1,jj)*omega(ii-1,jj);
end
if v(ii,jj)<0
laplvom(ii,jj) = v(ii,jj+1)*omega(ii,jj+1) - v(ii,jj)*omega(ii,jj);
else
laplvom(ii,jj) = v(ii,jj)*omega(ii,jj) -v(ii,jj-1)*omega(ii,jj-1);
end
end
end
omeganew = omega+dt.*( -lapluom/h -laplvom/h +nu.*laplomega); %omega for the next iteration
%making sure we begin with our circle again...
u(indcsToReset)=0; v(indcsToReset)=0;omeganew(indcsToReset)=0;
timestep = timestep+1;
% A = psinew(:);
treal=timestep*dt;
disp(treal)
end
for playin = 1:timestep
A = reshape(201,502,timestep);
M(playin) = A(:,:,playin);
end
movie(M)
figure(3); title('Steamline Plot with Karman Vortex Sheets');
contourf(psinew,100);set(gca,'DataAspectRatio',[1 1 1]);
toc
  5 Comments
Walter Roberson
Walter Roberson on 13 May 2014
I don't bother to run code to find all the errors. Tell us what error is occurring and on which line, and show us size() and class() of each variable mentioned on the line.
E
E on 13 May 2014
it occurs all in the while loop. all matrices are 201x502 double.
while error < max(max((F/4)*(psinew(1:r-2,2:c-2)+psi(3:r,2:c-2)+psi(2:r-1,3:c-1)+psinew(2:r-1,1:c-3)+(h^2).*omeganew(2:r-1,2:c-2)-4*psi(2:r-1,2:c-2)))) treal==0
psinew(2:r-1,2:c-2) = psi(2:r-1,2:c-2)+(F/4)*(psinew(1:r-2,2:c-2)+psi(3:r,2:c-2)+psi(2:r-1,3:c-1)+psinew(2:r-1,1:c-3)+(h^2).*omeganew(2:r-1,2:c-2)-4*psi(2:r-1,2:c-2));
psinew(:,c) = 2*psi(:,c-1)-psi(:,c-2);%--OUTFLOW BC
psinew(indcsToReset)=0;
psi=psinew;
%initialize vorticity at the wall
for ii = -100:100 %fluctuate i
for jj= 1:502 %fluctuate j
if (ii.^2 + (jj - j ./ 5).^2) < 625
omega(ii+101,jj)= -(2/(h^2))*(psinew(ii+1+101,jj)+psinew(ii-1+101,jj)+psinew(ii+101,jj+1)+psinew(ii+101,jj-1)-4*psinew(ii+101,jj));
end
end
end
%now find omeganew again
u(1:r-1,1:c-1) = (1/(2*h))*(psinew(1:r-1,2:c)- psinew(1:r-1,1:c-1)); u(:,c)=u(:,c-1); u(r,:)= u(r-1,:);
v(1:r-1,1:c-1) = (1/(2*h))*(psinew(1:r-1,1:c-1) - psinew(2:r,1:c-1));v(:,c)=v(:,c-1);v(r,:)=v(r-1,:);
laplomega=omeganew; omega = omeganew;
laplomega(2:r-1,2:c-2)= (1/(h^2))*(omega(1:r-2,2:c-2)+omega(3:r,2:c-2)+omega(2:r-1,3:c-1)+omega(2:r-1,1:c-3)-4*omega(2:r-1,2:c-2));
for ii = 2:r-1
for jj = 2:c-1
if u(ii,jj)<0
lapluom(ii,jj) = u(ii+1,j)*omega(ii+1,j)-u(ii,jj)*omega(ii,jj);
else
lapluom(ii,jj) = u(ii,jj)*omega(ii,jj)-u(ii-1,jj)*omega(ii-1,jj);
end
if v(ii,jj)<0
laplvom(ii,jj) = v(ii,jj+1)*omega(ii,jj+1) - v(ii,jj)*omega(ii,jj);
else
laplvom(ii,jj) = v(ii,jj)*omega(ii,jj) -v(ii,jj-1)*omega(ii,jj-1);
end
end
end
omeganew = omega+dt.*( -lapluom/h -laplvom/h +nu.*laplomega); %omega for the next iteration
%making sure we begin with our circle again...
u(indcsToReset)=0; v(indcsToReset)=0;omeganew(indcsToReset)=0; timestep = timestep+1; % A = psinew(:);
treal=timestep*dt; disp(treal) end

Sign in to comment.

Answers (0)

Categories

Find more on Thermal Analysis in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!