% This Matlab script solves the two-dimensional convection % equation using a finite volume algorithm on a mesh of (constant) % rectangular cells. The problem is assumed to be periodic and % have a constant velocity. close all; clear all; % Specify x range and number of points x0 = -2; x1 = 2; Nx = 40; % Specify y range and number of points y0 = -2; y1 = 2; Ny = 40; % 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)); [xmidg,ymidg] = ndgrid(xmid,ymid); % 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 final time tfinal = 10; % Set timestep CFL = 1.0; dt = CFL/(abs(u)/dx + abs(v)/dy); % Set initial condition to U0 = exp(-x^2 - 20*y^2) % Note: technically, we should average the initial % distribution in each cell but I chose to just set % the value of U in each control volume to the midpoint % value of U0. U = exp(-xmidg.^2 - 20*ymidg.^2); t = 0; % Loop until t > tfinal while (t < tfinal), % 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; % Copy U into Ubc Ubc( 1,2:Ny+1) = U(Nx, :); % Periodic bc Ubc(Nx+2,2:Ny+1) = U( 1, :); % Periodic bc Ubc(2:Nx+1, 1) = U( :,Ny); % Periodic bc Ubc(2:Nx+1,Ny+2) = U( :, 1); % Periodic bc % 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; % Forward Euler step U = U - (dt/A)*R; % Increment time t = t + dt; % Plot current solution Up = reshape(U,1,Nx*Ny); clf; [Hp] = patch(xp,yp,Up); set(Hp,'EdgeAlpha',0); axis('equal'); caxis([0,1]); colorbar; drawnow; end