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)
    load('cavity_variables.mat');

    % 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

    adu = 0.0;
    adv = 0.0;
    adt = 0.0;
    
    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.
            % Also adjust temperature for adiabatic boundary condition if
            % 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;

            % find the maximum changes for each equation - these are given by adu, adv, adt.

            ddu = abs(dpsi);
            ddv = abs(dw);
            ddt = abs(dtemp);
            
            if(ddu>adu)
                adu = ddu;
            end;
            if(ddv>adv)
                adv = ddv;
            end;
            if(ddt>adt)
                adt = ddt;
            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.
  
    if(adu<eps_psi && adv<eps_omega && adt<eps_temp)
        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.

Delete this entry (admin only).