# Demo entry 6742496

Cavity

Submitted by anonymous on May 18, 2018 at 05:49
Language: Matlab. Code size: 12.5 kB.

```function [x,y,re,pe,kount,nx,ny,relax1,relax2,relax3,a,psi,t,u,v,w,psm] = cavity(varargin)

% Program to solve for flow over a backward facing step.

%  clear all;
clc
% use this to restart from a previous solution

restart = 0;

% use this to plot streamfunction during iterations

iplot   = 0;

% counters ...

kount    = 0;

% set up the size of the grid
%
%   nx = cells in the x direction
%   ny = cells in the y direction

% define the number of cells in each direction
nx      =    64;
ny      =    64;

% define the Reynolds number and Peclet number here

re      = 100.0;
pe      = 100.0;

% define computational domain size

xlen     =  1.0;
ylen     =  1.0;

% define the relaxation parameters.
%
%   relax2 for streamfunction equation - this is well behaved so should
%   have relax2 > 1.
%   relax1 for the vorticity equation - this is nonlinear so require relax1
%   << 1 for reasonable Re.
%   relax3 for the temperature equation - again will require a small value
%   for this for reasonable Pe.

relax1   =  0.4;
relax2   =  1.5;
relax3   =  0.4;

% set convergence criteria for each equation.

eps_omega = 1.0e-6;
eps_psi   = 1.0e-6;
eps_temp  = 1.0e-6;

% set the parameters using from function argument list

for ii = 1:length(varargin)/2
if ~exist(varargin{2*ii-1})
error(['unknown parameter: ' varargin{2*ii-1}]);
end;
eval([varargin{2*ii-1} '=' num2str(varargin{2*ii}) ';']);
end;

fprintf('PROBLEM PARAMETERS FOR THIS RUN\n\n');
fprintf('  Cells in x direction        = %8d \n',nx);
fprintf('  Cells in y direction        = %8d \n\n',ny);

fprintf('  domain size in x            = %8.2f \n',xlen);
fprintf('  domain size in y            = %8.2f \n\n',ylen);

fprintf('  Reynolds number             = %8d \n',re);
fprintf('  Peclet number               = %8d \n\n',pe);

fprintf('  Omega relaxation parameter  = %8.2f \n',relax1);
fprintf('    Psi relaxation parameter  = %8.2f \n',relax2);
fprintf('   Temp relaxation parameter  = %8.2f \n\n',relax3);

fprintf('  Omega convergence criterion = %8.6f \n',eps_omega);
fprintf('    Psi convergence criterion = %8.6f \n',eps_psi);
fprintf('   Temp convergence criterion = %8.6f \n\n',eps_temp);

% Work out the constant cell size in x and y.

dx = xlen/nx;
dy = ylen/ny;

% and the viscosity and conduction coefficients: note length and velocity
% scales = 1

visc = 1/re;
cond = 1/pe;

% used later to simplify expressions.

anu  =        (2.0/dx/dx + 2.0/dy/dy);
rex  = visc * (2.0/dx/dx + 2.0/dy/dy);
ret  = cond * (2.0/dx/dx + 2.0/dy/dy);

% loop through setting up the initial fields with boundary condition
% correct if known

for i = 1:nx+1
for j = 1:ny+1
psi(i,j) = 0.0;
w(i,j)   = 0.0;
x(i,j)   = (i-1)*dx;
y(i,j)   = (j-1)*dy;
u(i,j)   = 0.0;
v(i,j)   = 0.0;
t(i,j)   = 0.0;
end;
% set boundary values for u at top and temperature at the bottom.
% INPUT REQUIRED HERE

u(i,ny+1)  = 1;
t(i,1   )  = 1;

end;

if(restart == 1)

% generate interpolation meshes to go from one grid to the other

n = size(psisave)
xinc = ((n(1)-1.0)/(nx))
yinc = ((n(2)-1.0)/(ny))

[xg,yg] = meshgrid(1:xinc:n(1),1:yinc:n(2));
[xf,yf] = meshgrid(1:n(1),1:n(2));
size(xg)

w    = (interp2(xf,yf,vortsave',xg,yg,'linear'))';
psi  = (interp2(xf,yf,psisave',xg,yg,'linear'))';
u    = (interp2(xf,yf,usave',xg,yg,'linear'))';
v    = (interp2(xf,yf,vsave',xg,yg,'linear'))';
t    = (interp2(xf,yf,tsave',xg,yg,'linear'))';

% make sure bcs are correct.

u(:,ny+1) = 1; u(:,1) = 0; u(nx+1,:) = 0.0; u(1,:) = 0;
psi(:,1) = 0; psi(:,ny+1) = 0; psi(1,:) = 0; psi(nx+1,:) = 0;
w(:,1) = 0;

if(iplot == 1)
psimax  = max(max(psi));
psimin  = min(min(psi));

for i = 1:10
vv(i) =  (i-1)*psimin/10.0;
end;
for i = 11:20
vv(i) = (i-10)*psimax/10;
end;

figure(1);
contourf(x,y,psi,vv);
daspect([1 1 1]);
xlim([0 1]);
xlabel('x');
ylabel('y');
lab = sprintf('Streamfunction \n Re = %6.1f,  Psi (min/max) = (%f, %f) \n', ...
re,psimin,psimax);
lab2 = sprintf('\n Grid = (%d, %d)',nx,ny);
title(strcat(lab,lab2));
drawnow;
end;

end;

% loop until convergence - limit total iterations to 100000 - this may not
% be enough?

for k = 1:1000000

psimin =  1000000.0;
psimax = -1000000.0;

% loop over all internal points.

for i = 2:nx
for j = 2:ny

% modify the neighbouring vorticity if at each boundary using
% internal values of the streamfunction and wall velocity.
% needed.

% left point
if(i ==  2)
w(i-1,j) = -2.0*((psi(i,j)-psi(i-1,j))+dx*v(i-1,j))/dx/dx;
t(i-1,j) = (4.0*t(i,j)-t(i+1,j))/3.0;
end;

% right point

% INPUT REQUIRED HERE .....

if(i == nx)
w(i+1,j) = -2.0*((psi(i,j)-psi(i+1,j))+dx*v(i+1,j))/dx/dx;
%t(i+1,j) = (4.0*t(i,j)-3*t(i-1,j));
t(i+1,j) = (4.0*t(i,j)-t(i-1,j))/3;
end;

% bottom point

if(j == 2)
w(i,j-1) = -2.0*((psi(i,j)-psi(i,j-1))-dy*u(i,j-1))/dy/dy;
end;

% top point

if(j == ny)
w(i,j+1) = -2.0*((psi(i,j)-psi(i,j+1))+dy*u(i,j+1))/dy/dy;
end;

% define the velocity given the streamfunction - updated during
% the loop

u(i,j) =  (psi(i  ,j+1)-psi(i  ,j-1))/2.0/dy;
v(i,j) = -(psi(i+1,j  )-psi(i-1,j  ))/2.0/dx;

% next Gauss-Seidel update - one point at a time

% INPUT REQUIRED HERE .....

wij   = -(2*visc*(1/dx/dx + 1/dy/dy))^-1*  (...
(u(i,j)/(2*dx))*(w(i+1,j)-w(i-1,j)) + ((v(i,j)/(2*dy))*(w(i,j+1)-w(i,j-1)))...
-visc*((1/dx/dx) * (w(i+1,j)+w(i-1,j)) + (1/dy/dy) * (w(i,j+1)+w(i,j-1))));

% determine the change from the previous value - required for
% relaxation

dw      = wij   - w(i,j);

%update the value.

w(i,j) = w(i,j) + relax1*dw;

% do the same for the streamfunction

psiij  = ((psi(i+1,j)+psi(i-1,j))/dx/dx + (psi(i,j+1)+psi(i,j-1))/dy/dy + w(i,j))/anu;

% determine the change from the previous value

dpsi   = psiij - psi(i,j);

psi(i,j) = psi(i,j) + relax2*dpsi;

% and the temperature field too.

% INPUT REQUIRED HERE .....

tij   = -(2*cond*(1/dx/dx + 1/dy/dy))^-1*  (...
(u(i,j)/(2*dx))*(t(i+1,j)-t(i-1,j)) + ((v(i,j)/(2*dy))*(t(i,j+1)-t(i,j-1)))...
-cond*(1/dx/dx * (t(i+1,j)+t(i-1,j)) + 1/dy/dy * (t(i,j+1)+t(i,j-1))));

% determine the change from the previous value

dtemp  = tij   - t(i,j);

%update the value.

t(i,j) = t(i,j) + relax3*dtemp;

% find maximum and minimum streamfunction over the domain and
% temperature where the streamfunction is minimum (tmin)

if(psi(i,j)>psimax)
psimax = psi(i,j);
end;
if(psi(i,j)<psimin)
psimin = psi(i,j);
tmin   = t(i,j);
end;

ddu = abs(dpsi);
ddv = abs(dw);
ddt = abs(dtemp);

end;
end;
end;

end;

end;

% print out the current state every hundred steps to check convergence.

psimax  = max(max(psi));
psimin  = min(min(psi));

kount       = kount + 1;
time(kount) = kount;
psm(kount)  = psimin;

if(mod(k,100)==0)

%       fprintf('k = %5d, max/min psi (%f, %f)\n',k, psimax, psimin);

for i = 1:10
vv(i) =  (i-1)*psimin/10.0;
end;
for i = 11:20
vv(i) = (i-10)*psimax/10;
end;

% plot streamlines as we go

if(iplot == 1)
figure(1);
contourf(x,y,psi,vv);
daspect([1 1 1]);
xlim([0 1]);
xlabel('x');
ylabel('y');
lab = sprintf('Streamfunction \n Re = %6.1f, Psi (min/max) = (%f, %f) \n', ...
re,psimin,psimax);
lab2 = sprintf('\n Grid = (%d, %d)',nx,ny);
title(strcat(lab,lab2));
drawnow;
end;

% save the convergence history

end;

if(mod(k,100)==0)
fprintf('k = %5d, max/min psi + tmin (%f, %f, %f), changes in psi/omega/temp (%f, %f, %f)\n',k, psimax, psimin, tmin, adu, adv, adt);
end;

% save the current solution for restart later...

if(mod(k,1000)==0)
system('copy cavity_variables.dat cavity_variables.mat');
vortsave = w;
psisave  = psi;
usave    = u;
vsave    = v;
tsave    = t;
save('cavity_variables.dat','vortsave','psisave','usave','vsave','tsave');
end;

if(mod(k,1000)==0 && iplot == 1)
figure(6);
plot(time,psm);
xlabel('iterations');
ylabel('Psimin');
title('convergence history');
end;

% break out of the loop if the maximum change in psi or omega or temperature is below
% some criterion.

break;
end;

end;

if(k <= 100)
fprintf('failed to converge\n');
a = 'not converged';
return;
end;

% set up a vector to contain values for plotting the streamfunction

for i = 1:10
vv(i) = -(i-1)*0.01;
end;
for i = 11:20
vv(i) = (i-10)*psimax/10;
end;

umax  = max(max(u));
umin  = min(min(u));
vmax  = max(max(v));
vmin  = min(min(v));
tempmax  = max(max(t));
tempmin  = min(min(t));
vortmax  = max(max(w));
vortmin  = min(min(w));

% plot streamlines
psimax  = max(max(psi));
psimin  = min(min(psi));

for i = 1:10
vv(i) =  (i-1)*psimin/10.0;
end;
for i = 11:20
vv(i) = (i-10)*psimax/10;
end;
figure(1);
contourf(x,y,psi,1000,'edgecolor','none');
daspect([1 1 1]);
xlim([0 1]);
xlabel('x');
ylabel('y');
lab = sprintf('Streamfunction \n Re = %6.1f,  Psi (min/max) = (%f, %f) \n', ...
re,psimin,psimax);
lab2 = sprintf('\n Grid = (%d, %d)',nx,ny);
title(strcat(lab,lab2));
colorbar;

% plot vorticity

figure(2);
%creating new w
wnew=w;
for i=length(w)
for j=length(w)
if w(i,j)<-2
wnew(i,j)=-2;
end
if w(i,j)>2
wnew(i,j)=2;
end
end
end
contourf(x,y,w,1000,'edgecolor','none');
daspect([1 1 1])
xlim([0 1]);
xlabel('x');
ylabel('y');
lab = sprintf('Vorticity \n Re = %6.1f,  Vorticity (min/max) = (%f, %f) \n', ...
re, vortmin,vortmax);
lab2 = sprintf('\n Grid = (%d, %d)',nx,ny);
title(strcat(lab,lab2));
colorbar;
caxis([-2, 2])

% plot temperature

figure(3);
contourf(x,y,t,1000,'edgecolor','none');
daspect([1 1 1])
xlim([0 1]);
xlabel('x');
ylabel('y');
lab = sprintf('Temperature \n Re = %6.1f,  Temperature (min/max) = (%f, %f) \n', ...
re, tempmin,tempmax);
lab2 = sprintf('\n Grid = (%d, %d)',nx,ny);
title(strcat(lab,lab2));
colorbar;

%plot velocity vectors

figure(4);
quiver(x,y,u,v);
daspect([1 1 1])
xlim([0 1]);
ylim([0 1]);
xlabel('x');
ylabel('y');
lab = sprintf('Velocity vectors \n Re = %6.1f,  U (min/max) = (%f, %f) \n', ...
re, umin,umax);
labb = sprintf('\n V = (%d, %d)',vmin,vmax);
lab2 = sprintf('\n Grid = (%d, %d)',nx,ny);
title(strcat(lab,labb,lab2));

figure(5);
plot(time,psm);
xlabel('iterations');
ylabel('Psimin');
title('convergence history');

% values we want
disp(min(min(abs(psi)>0)))
%find(psi=min(psi))

a = 'converged';
```

This snippet took 0.03 seconds to highlight.

Back to the Entry List or Home.