Solving Unsteady Navier stokes equations using Artificial compressibility method

clear all
close all
clc

%% Defining the problem domain
n_points = 101; % no_of_points
dom_length = 1;
h = dom_length/(n_points-1);
x = 0:h:dom_length; %X domain span
y = 0:h:dom_length; %Y domain span
dt = 0.05; %Time advancement
tau = 0.001; %Dual timestep
Re = 100; %Reynolds number
delta = 10; %Artificial compressibility
omega = pi/6;
time_val = 0;
final_time = 2*pi/omega;

%% Initializing the variables
%Final collocated variables
u_unsteady=0;
v_unsteady=0;
p_unsteady=1;

%Staggered variables
u(n_points+1,n_points)=0;
v(n_points,n_points+1)=0;
p(n_points+1,n_points+1)=1;
u(1,:)=2;

u_new(n_points+1,n_points)=0;
v_new(n_points,n_points+1)=0;
p_new(n_points+1,n_points+1)=1;
u_new(1,:)=2;

u_physical_old(n_points+1,n_points)=0;
v_physical_old(n_points,n_points+1)=0;
u_physical_old(1,:)=2;

%% Solving the governing equations
sol_index = 1;
iterations = 0;
error_req = 1e-6; %final required error residual

while time_val < final_time
    disp(['Solving for time = ' num2str(time_val) 's']) 
    error = 1;
    while error > 1e-6
    % x-momentum eq. - Interior
        for i = 2:n_points
            for j = 2:n_points - 1
                unsteady = (u(i,j) - u_physical_old(i,j))/dt;
                pressure = -(p(i,j+1) - p(i,j))/h;
                diffusion = (1/Re)*((u(i+1,j) - 2*u(i,j) + u(i-1,j))/(h*h) + (u(i,j+1) - 2*u(i,j) + u(i,j-1))/(h*h));
                advection_x = ((0.5*(u(i,j)+u(i,j+1)))^2 - (0.5*(u(i,j)+u(i,j-1)))^2)/h;
                advection_y = ((0.25*(u(i,j)+u(i-1,j))*(v(i-1,j)+v(i-1,j+1))) - (0.25*(u(i,j)+u(i+1,j))*(v(i,j)+v(i,j+1))))/h;
                u_new(i,j) = u(i,j) + tau*(-unsteady + diffusion - advection_x - advection_y + pressure);
            end
        end

        % x-momentum eq. - Boundary
        u_new(1,:) = 2*cos(omega*time_val) - u_new(2,:);
        u_new(n_points + 1,:) = -u_new(n_points,:);
        u_new(2:n_points,1) = 0;
        u_new(2:n_points,n_points) = 0;    

        % y-momentum eq. - Interior
        for i = 2:n_points - 1
            for j = 2:n_points
                unsteady = (v(i,j) - v_physical_old(i,j))/dt;
                pressure = -(p(i,j) - p(i+1,j))/h;
                diffusion = (1/Re)*((v(i+1,j) - 2*v(i,j) + v(i-1,j))/(h*h) + (v(i,j+1) - 2*v(i,j) + v(i,j-1))/(h*h));
                advection_y = ((0.5*(v(i,j)+v(i-1,j)))^2 - (0.5*(v(i,j)+v(i+1,j)))^2)/h;
                advection_x = ((0.25*(u(i,j)+u(i+1,j))*(v(i,j)+v(i,j+1))) - (0.25*(u(i,j-1)+u(i+1,j-1))*(v(i,j)+v(i,j-1))))/h;
                v_new(i,j) = v(i,j) + tau*(-unsteady + diffusion - advection_x - advection_y + pressure);
            end
        end

        % y-momentum eq. - Boundary
        v_new(:,1) = -v_new(:,2);
        v_new(:,n_points + 1) = -v_new(:,n_points);
        v_new(1,2:n_points) = 0;
        v_new(n_points,2:n_points) = 0;

        % Continuity eq. - Interior
        for i = 2:n_points
            for j = 2:n_points
                p_new(i,j) = p(i,j) - delta*tau*(u_new(i,j) - u_new(i,j-1) + v_new(i-1,j) - v_new(i,j))/h;
            end
        end

        % Continuity eq. - Boundary
        p_new(1,:) = p_new(2,:);
        p_new(n_points + 1,:) = p_new(n_points,:);
        p_new(:,1) = p_new(:,2);
        p_new(:,n_points + 1) = p_new(:,n_points);

        % Continuity residual as error measure
        error = 0;
        for i = 2:n_points - 1
            for j = 2:n_points - 1
                error = error + abs((u_new(i,j) - u_new(i,j-1) + v_new(i-1,j) - v_new(i,j))/h);
            end
        end

        % Finishing the iteration
        u = u_new;
        v = v_new;
        p = p_new;
        iterations = iterations + 1;
    end
    
    % Assigning the old physical time step to be newer one
    u_physical_old = u_new;
    v_physical_old = v_new;
    time_val = time_val + dt;
    
    % Saving the solution after every dt
    for i = 1:n_points
        for j = 1:n_points
            u_unsteady(sol_index,i,j) = 0.5*(u(i,j) + u(i+1,j));
            v_unsteady(sol_index,i,j) = 0.5*(v(i,j) + v(i,j+1));
            p_unsteady(sol_index,i,j) = 0.25*(p(i,j) + p(i,j+1) + p(i+1,j) + p(i+1,j+1));
        end
    end
    sol_index = sol_index + 1;
end

%% Unsteady contour plot
figure;
startx = 0.01:0.01:1;
starty = 0.5.*ones(length(startx),1);
count = 1;
for i = 20:20:180
    subplot(3,3,count)
    u_sample = u_unsteady(i,:,:);
    u_sample = reshape(u_sample, [size(u_sample,2) size(u_sample,3)]);
    x_dom = ((1:n_points)-1).*h;
    y_dom = 1-((1:n_points)-1).*h;
    [X,Y] = meshgrid(x_dom,y_dom);
    contourf(X,Y,u_sample, 21, 'LineStyle', 'none')
    caxis([-1 1])
    colorbar
    colormap('jet')
    xlabel('x')
    ylabel('y')
    count = count + 1;
    
%     streamline(X, Y, u_sample, v_sample, startx, starty)
end

Results: