% This Matlab script solves the two-dimensional convection % equation using a finite volume algorithm. % clear all; % Set flag for plotting doplot = 1; if (doplot == 1), figure(1); clf; figure(2); clf; figure(3); clf; end xymult = 2; % Specify x range and number of points x0 = -2; x1 = 2; Nx = xymult*20; % Specify y range and number of points y0 = -2; y1 = 2; Ny = xymult*20; % Construct mesh x = linspace(x0,x1,Nx+1); y = linspace(y0,y1,Ny+1); [xg,yg] = ndgrid(x,y); % Construct mesh needed for plotting xp = zeros(4,Nx*Ny); yp = zeros(4,Nx*Ny); n = 0; for j = 1:Ny, for i = 1:Nx, n = n + 1; xp(1,n) = x(i); yp(1,n) = y(j); xp(2,n) = x(i+1); yp(2,n) = y(j); xp(3,n) = x(i+1); yp(3,n) = y(j+1); xp(4,n) = x(i); yp(4,n) = y(j+1); end end % Calculate midpoint values in each control volume xmid = 0.5*(x(1:Nx) + x(2:Nx+1)); ymid = 0.5*(y(1:Ny) + y(2:Ny+1)); % Calculate cell size in control volumes (assumed equal) dx = x(2) - x(1); dy = y(2) - y(1); A = dx*dy; % Set velocity u = 1; v = -1; % Set tolerances on iteration and residual tolerance ntol = 3000; Rtol = -13; % Set timestep CFL = 1.0; dt = CFL/(abs(u)/dx + abs(v)/dy); % Set initial conditions U = zeros(Nx,Ny); n = 0; Rmax = 0; % Loop until one of the termination criterion are met while ((n < ntol) & (Rmax > Rtol)), % The following implement the bc's by creating a larger array % for U and putting the appropriate values in the first and last % columns or rows to set the correct bc's Ubc(2:Nx+1,2:Ny+1) = U; Ubc( 1,2:Ny+1) = exp(-10*(ymid-1).^2); % Ubc(Nx+2,2:Ny+1) = zeros(1, Ny); % Ubc(2:Nx+1, 1) = zeros(Nx, 1); % Ubc(2:Nx+1,Ny+2) = exp(-10*(y1-1)^2); % % Calculate the flux at each interface % First the i interfaces F = 0.5* u *( Ubc(2:Nx+2,2:Ny+1) + Ubc(1:Nx+1,2:Ny+1)) ... - 0.5*abs(u)*( Ubc(2:Nx+2,2:Ny+1) - Ubc(1:Nx+1,2:Ny+1)); % Now the j interfaces G = 0.5* v *( Ubc(2:Nx+1,2:Ny+2) + Ubc(2:Nx+1,1:Ny+1)) ... - 0.5*abs(v)*( Ubc(2:Nx+1,2:Ny+2) - Ubc(2:Nx+1,1:Ny+1)); % Add contributions to residuals from fluxes R = (F(2:Nx+1,:) - F(1:Nx,:))*dy + (G(:,2:Ny+1) - G(:,1:Ny))*dx; % Calculate maximum residual Rmax = log10(max(max(abs(R)))); % Forward Euler step U = U - (dt/A)*R; % Increment iteration n = n + 1; fprintf('Iter: %i, Log10 Max Res = %f\n',n,Rmax); if (doplot == 1), % Plot current solution Up = reshape(U,1,Nx*Ny); figure(1); clf; [Hp] = patch(xp,yp,Up); set(Hp,'EdgeAlpha',0); axis('equal'); title(sprintf('Log10 Max Res = %f\n',Rmax)); xlabel('x'); ylabel('y'); caxis([0,1]); colorbar; % Plot max residual figure(2); hold on; plot(n,Rmax,'*'); hold off; xlabel('Iteration'); ylabel('Log10 Rmax'); % Plot solution at bottom of domain figure(3); stairs(x,[U(:,1)', U(Nx,1)]); xlabel('x'); ylabel('U'); axis([x0,x1,0,1]); title(sprintf('U at y = %f\n',y0)); drawnow; end end