FINAL INDEPENDENT PROJECT
To put it simply as throughout the corse there were no switche case programs I wanted in incoperate switch case function in my program
AIM:- To find the solutions for various items using the navier stokes equation there are 11 cases programed each with varying levels of complecity.
Note:- for the purposes of avoiding too many errors and obtain the right outputs. fixed variables are used for all the cases. however we can change the input parameters.
clear
close all
clc
%% Input variables:
fprintf('There are 10 cases to chose from \n\t case 1 \n One dimensional convection diffusion \n\t case 2 \n second derivative \n\t case 3 \n flow through channel \n\t case 4 \n conventional diffusion upwind \n\t case 5 \n conventional diffusion \n\t case 6 \n two dimensional diffusion \n\t case 7 \n two dimensional LDC ACM \n\t case 8 \n two dimensional LDC simple alternate \n\t case 9 \n two dimensional LDC \n\t case 10 \n two dimensional UNS diffusion plot \n\t case 11 \n two dimensional unsteady LDC ACM \n\n\n');
Navierstroke = input('Enter the Problem number needed to be solved\n');
%% Using A Switch case
switch Navierstroke
%% CASE 1:
case 1
%% Defining the problem
n_points = 101;
dom_length = 1;
h = dom_length/(n_points-1);
x = 0:h:dom_length;
phi(1) = 10;
phi(n_points) = 20;
phi_new(1) = 10;
phi_new(n_points) = 20;
rho = 1;
u = 1;
gamma = 1;
P = rho*u*dom_length/gamma;
%% Exact solution
phi_exact = phi(1) + (phi(n_points)-phi(1)).*(exp(P.*x/dom_length) - 1)./(exp(P) - 1);
figure(1);
plot(x, phi_exact, 'rp')
%% Differencing scheme
error = 1;
iterations = 0;
while error > 1e-7
for i = 2:n_points-1
a_E = gamma/h - rho*u/2;
a_W = gamma/h + rho*u/2;
a_P = rho*u/2 - rho*u/2 + gamma/h + gamma/h;
phi_new(i) = (a_E*phi(i+1) + a_W*phi(i-1))/a_P;
end
iterations = iterations + 1;
error = 0;
for i = 2:n_points-1
error = error + abs(phi(i) - phi_new(i));
end
phi = phi_new;
end
figure(1);
hold on
plot(x, phi, 'b', 'LineWidth', 1)
xlabel('x')
ylabel('Phi')
legend('Exact solution', 'Numerical scheme')
%% CASE 2:
case 2
%% Defining the mesh
n_points = 5;
dom_size = 1;
h = dom_size/(n_points - 1);
%% Initialising the problem
y(1) = 0;
y(n_points) = 1;
y_new(1) = 0;
y_new(n_points) = 1;
error_mag = 1;
error_req = 1e-7;
iterations = 0;
%% Calculations
while error_mag > error_req
for i = 2:(n_points-1)
y_new(i) = 0.5.*(y(i-1) + y(i+1));
iterations = iterations + 1;
end
% Calculation of error magnitude
error_mag = 0;
for i = 2:(n_points-1)
error_mag = error_mag + abs(y(i) - y_new(i));
end
% Assigning new to be old
y = y_new;
end
%% Plotting
x_dom = ((1:n_points)-1).*h;
figure; plot(x_dom, y)
%% CASE 3:
case 3
%% Defining the problem domain
n_points_y = 11;
n_points_x = 101;
dom_length_x = 10;
dom_length_y = 1;
h = dom_length_y/(n_points_y-1);
x = 0:h:dom_length_x;
y = 0:h:dom_length_y;
mu = 0.01;
rho = 1;
% Under-relaxation factors
alpha = 0.1;
alpha_p = 0.8;
%% Initializing the variables
% Final collocated variables
u_final(n_points_y,n_points_x) = 0;
v_final(n_points_y,n_points_x) = 0;
p_final(n_points_y,n_points_x) = 0;
u_final(:,1) = 1;
u_final(:,n_points_x) = 1;
% Staggered variables
u(n_points_y+1,n_points_x) = 0;
ustar(n_points_y+1,n_points_x) = 0;
d_u(n_points_y+1,n_points_x) = 0;
v(n_points_y,n_points_x+1) = 0;
vstar(n_points_y,n_points_x+1) = 0;
d_v(n_points_y,n_points_x+1) = 0;
p(n_points_y+1,n_points_x+1) = 0;
pstar(n_points_y+1,n_points_x+1) = 0;
pcor(n_points_y+1,n_points_x+1) = 0;
b(n_points_y+1,n_points_x+1) = 0;
u(:,1) = 1;
u(:,n_points_x) = 1;
ustar(:,1) = 1;
ustar(:,n_points_x) = 1;
u_new(n_points_y+1,n_points_x) = 0;
v_new(n_points_y,n_points_x+1) = 0;
p_new(n_points_y+1,n_points_x+1) = 0;
u_new(:,1) = 1;
u_new(:,n_points_x) = 1;
%% Solving the governing equations
error = 1;
iterations = 0;
error_req = 1e-7;
figure(1);
while error > error_req
% x-momentum eq. interior
for i = 2:n_points_y
for j = 2:n_points_x-1
a_E = (mu/h) - (rho*(u(i,j)+u(i,j+1))/4);
a_W = (mu/h) + (rho*(u(i,j)+u(i,j-1))/4);
a_N = (mu/h) - (rho*(v(i-1,j)+v(i-1,j+1))/4);
a_S = (mu/h) + (rho*(v(i,j)+v(i,j+1))/4);
a_P = a_E + a_W + a_N + a_S + rho*((u(i,j)+u(i,j+1))/2 - (u(i,j)+u(i,j-1))/2 + (v(i-1,j)+v(i-1,j+1))/2 - (v(i,j)+v(i,j+1))/2);
d_u(i,j) = h/a_P;
ustar_mid = (a_E*u(i,j+1) + a_W*u(i,j-1) + a_N*u(i-1,j) + a_S*u(i+1,j) + (p(i,j)-p(i,j+1))*h)/a_P;
ustar(i,j) = (1-alpha)*u(i,j) + alpha*ustar_mid;
end
end
% x-momentum eq. boundary
ustar(:,1) = ustar(:,n_points_x-1);
ustar(:,n_points_x) = ustar(:,2);
ustar(1,2:n_points_x) = -ustar(2,2:n_points_x);
ustar(n_points_y+1,2:n_points_x) = -ustar(n_points_y,2:n_points_x);
%y-momentum eq. interior
for i = 2:n_points_y-1
for j = 2:n_points_x
a_E = (mu/h) - (rho*(u(i,j)+u(i+1,j))/4);
a_W = (mu/h) + (rho*(u(i,j-1)+u(i+1,j-1))/4);
a_N = (mu/h) - (rho*(v(i,j)+v(i-1,j))/4);
a_S = (mu/h) + (rho*(v(i,j)+v(i+1,j))/4);
a_P = a_E + a_W + a_N + a_S + rho*((u(i,j)+u(i+1,j))/2 - (u(i,j-1)+u(i+1,j-1))/2 + (v(i,j)+v(i-1,j))/2 - (v(i,j)+v(i+1,j))/2);
d_v(i,j) = h/a_P;
vstar_mid = (a_E*v(i,j+1) + a_W*v(i,j-1) + a_N*v(i-1,j) + a_S*v(i+1,j) + (p(i+1,j)-p(i,j))*h)/a_P;
vstar(i,j) = (1-alpha)*v(i,j) + alpha*vstar_mid;
end
end
% y-momentum eq. boundary
vstar(1,2:n_points_x) = 0;
vstar(n_points_y,2:n_points_x) = 0;
vstar(:,1) = vstar(:,n_points_x);
vstar(:,n_points_x+1) = vstar(:,2);
% Zeroing the corrections to begin with
pcor(1:n_points_y+1,1:n_points_x+1) = 0;
% Continuity eq. aka pressure correction - interior
for i = 2:n_points_y
for j = 2:n_points_x
a_E = rho*h*d_u(i,j);
a_W = rho*h*d_u(i,j-1);
a_N = rho*h*d_v(i-1,j);
a_S = rho*h*d_v(i,j);
a_P = a_E + a_W + a_N + a_S;
b(i,j) = rho*(h*(ustar(i,j-1)-ustar(i,j)) + h*(vstar(i,j)-vstar(i-1,j)));
pcor(i,j) = (a_E*pcor(i,j+1) + a_W*pcor(i,j-1) + a_N*pcor(i-1,j) + a_S*pcor(i+1,j) + b(i,j))/a_P;
end
end
% Correcting the pressure field
for i = 2:n_points_y
for j = 2:n_points_x
p_new(i,j) = p(i,j) + alpha_p*pcor(i,j);
end
end
% Continuity eq boundary
p_new(1:n_points_y,n_points_x+1) = p_new(1:n_points_y,2);
p_new(1:n_points_y,1) = p_new(1:n_points_y,n_points_x);
p_new(1,1:n_points_x) = p_new(2,1:n_points_x);
p_new(n_points_y+1,1:n_points_x) = p_new(n_points_y,1:n_points_x);
% Correcting the velocities
for i = 2:n_points_y
for j = 2:n_points_x-1
u_new(i,j) = ustar(i,j) + d_u(i,j)*(pcor(i,j)-pcor(i,j+1));
end
end
% x-momentum eq. boundary
u_new(:,1) = u_new(:,n_points_x-1);
u_new(:,n_points_x) = u_new(:,2);
u_new(1,2:n_points_x) = -u_new(2,2:n_points_x);
u_new(n_points_y+1,2:n_points_x) = -u_new(n_points_y,2:n_points_x);
for i = 2:n_points_y-1
for j = 2:n_points_x
v_new(i,j) = vstar(i,j) + d_v(i,j)*(pcor(i+1,j)-pcor(i,j));
end
end
% y-momentum eq. boundary
v_new(1,:) = 0;
v_new(n_points_y,:) = 0;
v_new(2:n_points_y-1,1) = v_new(2:n_points_y-1,n_points_x);
v_new(2:n_points_y-1,n_points_x+1) = v_new(2:n_points_y-1,2);
% continuity eq. as error measure
error = 0;
for i = 2:n_points_y
for j = 2:n_points_x
error = error + abs(b(i,j));
end
end
% error monitoring after every few timesteps
if(rem(iterations,1000)) == 0
figure(1);
semilogy(iterations, error, '-ko')
hold on
xlabel('Iterations');
ylabel('Residual Error');
end
% Finishing the iteration
u = u_new;
v = v_new;
p = p_new;
iterations = iterations + 1;
end
% Mapping the staggered variables to collocated variables
for i = 1:n_points_y
for j = 1:n_points_x
u_final(i,j) = 0.5*(u(i,j) + u(i+1,j));
v_final(i,j) = 0.5*(v(i,j) + v(i,j+1));
p_final(i,j) = 0.25*(p(i,j) + p(i,j+1) + p(i+1,j) + p(i+1,j+1));
end
end
%% Contour and vector visuals.
x_dom = ((1:n_points_x)-1).*h;
y_dom = 1-((1:n_points_y)-1).*h;
[X,Y] = meshgrid(x_dom,y_dom);
figure(21);
contourf(X,Y, u_final, 21, 'LineStyle', 'none')
colorbar
colormap('jet')
xlabel('x')
ylabel('y')
title('U- Velocity')
daspect([1 0.5 2])
figure(22);
contourf(X,Y, p_final, 21, 'LineStyle', 'none')
colorbar
colormap('jet')
xlabel('x')
ylabel('y')
title('Pressure')
daspect([1 0.5 2])
figure(23);
hold on
quiver(X, Y, u_final, v_final, 1, 'k')
daspect([1 0.5 2])
%% CASE 4:
case 4
%% 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
%% Initializing the problem
T(1,1:n_points) = 1;
T(1:n_points,1) = 1;
T_new(1,1:n_points) = 1;
T_new(1:n_points,1) = 1;
rho = 1;
u = 1;
v = 1;
gamma = 0.001;
P = rho*u*h/gamma;
%% Upwind scheme
error = 1;
iterations = 0;
while error > 1e-7
for i = 2:n_points-1
for j = 2:n_points-1
a_E = gamma;
a_W = gamma + rho*u*h;
a_N = gamma;
a_S = gamma + rho*v*h;
a_P = rho*u*h + rho*v*h + gamma + gamma + gamma + gamma;
T_new(i,j) = (a_E*T(i+1,j) + a_W*T(i-1,j) + a_N*T(i,j-1) + a_S*T(i,j+1))/a_P;
end
end
iterations = iterations + 1;
error = 0;
for i = 2:n_points-1
for j = 2:n_points-1
error = error + abs(T(i,j) - T_new(i,j));
end
end
T = T_new;
end
%% Plotting
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,T, 50)
colorbar
%% Centerline temperature
figure;
plot(1-y,T(:,(n_points+1)/2),'--o')
%% CASE 5:
case 5
%% 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
%% Initializing the problem
T(1,1:n_points) = 1;
T(1:n_points,1) = 1;
T_new(1,1:n_points) = 1;
T_new(1:n_points,1) = 1;
rho = 1;
u = 1;
v = 1;
gamma = 0.0025;
P = rho*u*h/gamma
%% Differencing scheme
error = 1;
iterations = 0;
while error > 1e-7
for i = 2:n_points-1
for j = 2:n_points-1
a_E = gamma - rho*u*h/2;
a_W = gamma + rho*u*h/2;
a_N = gamma - rho*v*h/2;
a_S = gamma + rho*v*h/2;
a_P = rho*u*h/2 - rho*u*h/2 + rho*v*h/2 - rho*v*h/2 + gamma + gamma + gamma + gamma;
T_new(i,j) = (a_E*T(i+1,j) + a_W*T(i-1,j) + a_N*T(i,j-1) + a_S*T(i,j+1))/a_P;
end
end
iterations = iterations + 1;
error = 0;
for i = 2:n_points-1
for j = 2:n_points-1
error = error + abs(T(i,j) - T_new(i,j));
end
end
T = T_new;
end
%% Plotting
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,T, 50)
colorbar
%% Centerline temperature
figure;
plot(1-y,T(:,(n_points+1)/2),'--o')
%% CASE 6:
case 6
%% Defining the mesh
n_points = 101;
dom_size = 1;
h = dom_size/(n_points - 1);
dt = 0.000001;
alpha = dt/(h*h);
%% Initialising the problem
y(n_points, n_points) = 0;
y(1,:) = 1;
y_new(n_points, n_points) = 0;
y_new(1,:) = 1;
error_mag = 1;
error_req = 1e-6;
iterations = 0;
%Tracking the error magnitude
error_track = 0;
%% Calculations
while error_mag > error_req
for i = 2:(n_points-1)
for j = 2:(n_points-1)
y_new(i,j) = y(i,j) + alpha.*(y(i-1,j) + y(i+1,j) + y(i,j-1) + y(i,j+1) - 4*y(i,j));
end
end
iterations = iterations + 1;
% Calculation of error magnitude
error_mag = 0;
for i = 2:(n_points-1)
for j = 2:(n_points-1)
error_mag = error_mag + abs(y(i,j) - y_new(i,j));
error_track(iterations) = error_mag;
end
end
if rem(iterations, 1000) == 0
iterations
error_mag
end
% Assigning new to be old
y = y_new;
end
%% Plotting
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,y, 12)
colorbar
% Plot the error with time
figure;
time = dt.*(1:iterations);
plot(time, error_track)
%% Concept of subplots
figure;
subplot(2,1,1)
plot(time, error_track)
subplot(2,1,2)
contourf(X,Y,y, 12)
colorbar
%% CASE 7:
case 7
%% Defining the problem domain
n_points = 51; % 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.001; %Time advancement
Re = 100; %Reynolds number
delta = 4.5; %Artificial compressibility
%% Initializing the variables
%Final collocated variables
u_final(n_points,n_points)=0;
v_final(n_points,n_points)=0;
p_final(n_points,n_points)=1;
u_final(1,:) = 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;
%% Solving the governing equations
error = 1;
iterations = 0;
error_req = 1e-6; %final required error residual
figure(1); %for error monitoring
while error > error_req
% x-momentum eq. - Interior
for i = 2:n_points
for j = 2:n_points - 1
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) + dt*(diffusion - advection_x - advection_y + pressure);
end
end
% x-momentum eq. - Boundary
u_new(1,:) = 2 - 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
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) + dt*(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*dt*(u(i,j) - u(i,j-1) + v(i-1,j) - v(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
% Error monitoring after every few timesteps
if(rem(iterations, 5000)) == 0
figure(1);
semilogy(iterations, error, '-ko')
hold on
xlabel('Iterations')
ylabel('Residual Error')
end
% Finishing the iteration
u = u_new;
v = v_new;
p = p_new;
iterations = iterations + 1;
end
% After the converged solution, we map the staggered variables to
% collocated variables
for i = 1:n_points
for j = 1:n_points
u_final(i,j) = 0.5*(u(i,j) + u(i+1,j));
v_final(i,j) = 0.5*(v(i,j) + v(i,j+1));
p_final(i,j) = 0.25*(p(i,j) + p(i,j+1) + p(i+1,j) + p(i+1,j+1));
end
end
%% Centreline u validation - Comparison with benchmark
figure(11);
plot(u_final(:,(n_points+1)/2),1-y, 'LineWidth', 1)
tab_data = xlsread('plot_u_y_Ghia100.csv', 'A2:B18');
y_ghia = tab_data(:,1);
u_ghia = tab_data(:,2);
figure(11); hold on
plot(u_ghia, y_ghia, 'o', 'LineWidth', 1)
xlabel('u')
ylabel('y')
legend('Numerical', 'Benchmark', 'location', 'southeast')
%% Contour and vector visuals.
x_dom = ((1:n_points)-1).*h;
y_dom = 1-((1:n_points)-1).*h;
[X,Y] = meshgrid(x_dom,y_dom);
figure(21);
contourf(X,Y,v_final, 21, 'LineStyle', 'none')
colorbar
colormap('jet')
xlabel('x')
ylabel('y')
figure(22);
hold on
quiver(X, Y, u_final, v_final, 5, 'k')
%% CASE 8:
case 8
%% Defining the problem domain
n_points = 51; % 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
Re = 1000; %Reynolds number
nu = 1/Re;
% Under-relaxation factors
alpha = 0.01;
alpha_p = 0.8;
%% Initializing the variables
%Final collocated variables
u_final(n_points,n_points)=0;
v_final(n_points,n_points)=0;
p_final(n_points,n_points)=1;
u_final(1,:) = 1;
%Staggered variables
u(n_points+1,n_points)=0;
u_star(n_points+1,n_points)=0;
d_e(n_points+1,n_points)=0;
v(n_points,n_points+1)=0;
v_star(n_points,n_points+1)=0;
d_n(n_points,n_points+1)=0;
p(n_points+1,n_points+1)=1;
p_star(n_points+1,n_points+1)=1;
pc(n_points+1,n_points+1)=0;
b(n_points+1,n_points+1)=0;
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;
%% Solving the governing equations
error = 1;
iterations = 0;
error_req = 1e-7; %final required error residual
figure(1); %for error monitoring
while error > error_req
% x-momentum eq. - Interior
for i = 2:n_points
for j = 2:n_points - 1
u_E = 0.5*(u(i,j) + u(i,j+1));
u_W = 0.5*(u(i,j) + u(i,j-1));
v_N = 0.5*(v(i-1,j) + v(i-1,j+1));
v_S = 0.5*(v(i,j) + v(i,j+1));
a_E = -0.5*u_E*h + nu;
a_W = 0.5*u_W*h + nu;
a_N = -0.5*v_N*h + nu;
a_S = 0.5*v_S*h + nu;
a_e = 0.5*u_E*h - 0.5*u_W*h + 0.5*v_N*h - 0.5*v_S*h + 4*nu;
A_e = -h;
d_e(i,j) = A_e/a_e;
u_star_mid = (a_E*u(i,j+1) + a_W*u(i,j-1) + a_N*u(i-1,j) + a_S*u(i+1,j))/a_e + d_e(i,j)*(p(i,j+1) - p(i,j));
u_star(i,j) = (1-alpha)*u(i,j) + alpha*u_star_mid;
end
end
% x-momentum eq. - Boundary
u_star(1,:) = 2 - u_star(2,:);
u_star(n_points + 1,:) = -u_star(n_points,:);
u_star(2:n_points,1) = 0;
u_star(2:n_points,n_points) = 0;
% y-momentum eq. - Interior
for i = 2:n_points - 1
for j = 2:n_points
u_E = 0.5*(u(i,j) + u(i+1,j));
u_W = 0.5*(u(i,j-1) + u(i+1,j-1));
v_N = 0.5*(v(i-1,j) + v(i,j));
v_S = 0.5*(v(i,j) + v(i+1,j));
a_E = -0.5*u_E*h + nu;
a_W = 0.5*u_W*h + nu;
a_N = -0.5*v_N*h + nu;
a_S = 0.5*v_S*h + nu;
a_n = 0.5*u_E*h - 0.5*u_W*h + 0.5*v_N*h - 0.5*v_S*h + 4*nu;
A_n = -h;
d_n(i,j) = A_n/a_n;
v_star_mid = (a_E*v(i,j+1) + a_W*v(i,j-1) + a_N*v(i-1,j) + a_S*v(i+1,j))/a_n + d_n(i,j)*(p(i,j) - p(i+1,j));
v_star(i,j) = (1-alpha)*v(i,j) + alpha*v_star_mid;
end
end
% y-momentum eq. - Boundary
v_star(:,1) = -v_star(:,2);
v_star(:,n_points + 1) = -v_star(:,n_points);
v_star(1,2:n_points) = 0;
v_star(n_points,2:n_points) = 0;
% Zeroing the corrections to begin with
pc(1:n_points+1,1:n_points+1)=0;
% Continuity equation a.k.a. pressure correction - Interior
for i = 2:n_points
for j = 2:n_points
a_E = -d_e(i,j)*h;
a_W = -d_e(i,j-1)*h;
a_N = -d_n(i-1,j)*h;
a_S = -d_n(i,j)*h;
a_P = a_E + a_W + a_N + a_S;
b(i,j) = -(u_star(i,j) - u_star(i,j-1))*h + (v_star(i,j) - v_star(i-1,j))*h;
pc(i,j) = (a_E*pc(i,j+1) + a_W*pc(i,j-1) + a_N*pc(i-1,j) + a_S*pc(i+1,j) + b(i,j))/a_P;
end
end
% Correcting the pressure field
for i = 2:n_points
for j = 2:n_points
p_new(i,j) = p(i,j) + alpha_p*pc(i,j);
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);
% Correcting the velocities
for i = 2:n_points
for j = 2:n_points - 1
u_new(i,j) = u_star(i,j) + d_e(i,j)*(pc(i,j+1) - pc(i,j));
end
end
% x-momentum eq. - Boundary
u_new(1,:) = 2 - 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;
for i = 2:n_points - 1
for j = 2:n_points
v_new(i,j) = v_star(i,j) + d_n(i,j)*(pc(i,j) - pc(i+1,j));
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 residual as error measure
error = 0;
for i = 2:n_points
for j = 2:n_points
error = error + abs(b(i,j));
end
end
% Error monitoring after every few timesteps
if(rem(iterations, 1000)) == 0
figure(1);
semilogy(iterations, error, '-ko')
hold on
xlabel('Iterations')
ylabel('Residual Error')
end
% Finishing the iteration
u = u_new;
v = v_new;
p = p_new;
iterations = iterations + 1;
end
% After the converged solution, we map the staggered variables to
% collocated variables
for i = 1:n_points
for j = 1:n_points
u_final(i,j) = 0.5*(u(i,j) + u(i+1,j));
v_final(i,j) = 0.5*(v(i,j) + v(i,j+1));
p_final(i,j) = 0.25*(p(i,j) + p(i,j+1) + p(i+1,j) + p(i+1,j+1));
end
end
%% Centreline u validation - Comparison with benchmark
figure(11);
plot(u_final(:,(n_points+1)/2),1-y, 'LineWidth', 1)
tab_data = xlsread('plot_u_y_Ghia100.csv', 'A2:B18');
y_ghia = tab_data(:,1);
u_ghia = tab_data(:,2);
figure(11); hold on
plot(u_ghia, y_ghia, 'o', 'LineWidth', 1)
xlabel('u')
ylabel('y')
legend('Numerical', 'Benchmark', 'location', 'southeast')
%% Contour and vector visuals.
x_dom = ((1:n_points)-1).*h;
y_dom = 1-((1:n_points)-1).*h;
[X,Y] = meshgrid(x_dom,y_dom);
figure(21);
contourf(X,Y,u_final, 21, 'LineStyle', 'none')
colorbar
colormap('jet')
xlabel('x')
ylabel('y')
figure(22);
hold on
quiver(X, Y, u_final, v_final, 5, 'k')
%% CASE 9:
case 9
%% Defining the problem domain
n_points = 51; % 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
Re = 100; %Reynolds number
nu = 1/Re;
% Under-relaxation factors
alpha = 0.8;
alpha_p = 0.8;
%% Initializing the variables
%Final collocated variables
u_final(n_points,n_points)=0;
v_final(n_points,n_points)=0;
p_final(n_points,n_points)=1;
u_final(1,:) = 1;
%Staggered variables
u(n_points+1,n_points)=0;
u_star(n_points+1,n_points)=0;
d_e(n_points+1,n_points)=0;
v(n_points,n_points+1)=0;
v_star(n_points,n_points+1)=0;
d_n(n_points,n_points+1)=0;
p(n_points+1,n_points+1)=1;
p_star(n_points+1,n_points+1)=1;
pc(n_points+1,n_points+1)=0;
b(n_points+1,n_points+1)=0;
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;
%% Solving the governing equations
error = 1;
iterations = 0;
error_req = 1e-7; %final required error residual
figure(1); %for error monitoring
while error > error_req
% x-momentum eq. - Interior
for i = 2:n_points
for j = 2:n_points - 1
u_E = 0.5*(u(i,j) + u(i,j+1));
u_W = 0.5*(u(i,j) + u(i,j-1));
v_N = 0.5*(v(i-1,j) + v(i-1,j+1));
v_S = 0.5*(v(i,j) + v(i,j+1));
a_E = -0.5*u_E*h + nu;
a_W = 0.5*u_W*h + nu;
a_N = -0.5*v_N*h + nu;
a_S = 0.5*v_S*h + nu;
a_e = 0.5*u_E*h - 0.5*u_W*h + 0.5*v_N*h - 0.5*v_S*h + 4*nu;
A_e = -h;
d_e(i,j) = A_e/a_e;
u_star(i,j) = (a_E*u(i,j+1) + a_W*u(i,j-1) + a_N*u(i-1,j) + a_S*u(i+1,j))/a_e + d_e(i,j)*(p(i,j+1) - p(i,j));
end
end
% x-momentum eq. - Boundary
u_star(1,:) = 2 - u_star(2,:);
u_star(n_points + 1,:) = -u_star(n_points,:);
u_star(2:n_points,1) = 0;
u_star(2:n_points,n_points) = 0;
% y-momentum eq. - Interior
for i = 2:n_points - 1
for j = 2:n_points
u_E = 0.5*(u(i,j) + u(i+1,j));
u_W = 0.5*(u(i,j-1) + u(i+1,j-1));
v_N = 0.5*(v(i-1,j) + v(i,j));
v_S = 0.5*(v(i,j) + v(i+1,j));
a_E = -0.5*u_E*h + nu;
a_W = 0.5*u_W*h + nu;
a_N = -0.5*v_N*h + nu;
a_S = 0.5*v_S*h + nu;
a_n = 0.5*u_E*h - 0.5*u_W*h + 0.5*v_N*h - 0.5*v_S*h + 4*nu;
A_n = -h;
d_n(i,j) = A_n/a_n;
v_star(i,j) = (a_E*v(i,j+1) + a_W*v(i,j-1) + a_N*v(i-1,j) + a_S*v(i+1,j))/a_n + d_n(i,j)*(p(i,j) - p(i+1,j));
end
end
% y-momentum eq. - Boundary
v_star(:,1) = -v_star(:,2);
v_star(:,n_points + 1) = -v_star(:,n_points);
v_star(1,2:n_points) = 0;
v_star(n_points,2:n_points) = 0;
% Zeroing the corrections to begin with
pc(1:n_points+1,1:n_points+1)=0;
% Continuity equation a.k.a. pressure correction - Interior
for i = 2:n_points
for j = 2:n_points
a_E = -d_e(i,j)*h;
a_W = -d_e(i,j-1)*h;
a_N = -d_n(i-1,j)*h;
a_S = -d_n(i,j)*h;
a_P = a_E + a_W + a_N + a_S;
b(i,j) = -(u_star(i,j) - u_star(i,j-1))*h + (v_star(i,j) - v_star(i-1,j))*h;
pc(i,j) = (a_E*pc(i,j+1) + a_W*pc(i,j-1) + a_N*pc(i-1,j) + a_S*pc(i+1,j) + b(i,j))/a_P;
end
end
% Correcting the pressure field
for i = 2:n_points
for j = 2:n_points
p_new(i,j) = p(i,j) + alpha_p*pc(i,j);
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);
% Correcting the velocities
for i = 2:n_points
for j = 2:n_points - 1
u_new(i,j) = u_star(i,j) + alpha*d_e(i,j)*(pc(i,j+1) - pc(i,j));
end
end
% x-momentum eq. - Boundary
u_new(1,:) = 2 - 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;
for i = 2:n_points - 1
for j = 2:n_points
v_new(i,j) = v_star(i,j) + alpha*d_n(i,j)*(pc(i,j) - pc(i+1,j));
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 residual as error measure
error = 0;
for i = 2:n_points
for j = 2:n_points
error = error + abs(b(i,j));
end
end
% Error monitoring after every few timesteps
if(rem(iterations, 1000)) == 0
figure(1);
semilogy(iterations, error, '-ko')
hold on
xlabel('Iterations')
ylabel('Residual Error')
end
% Finishing the iteration
u = u_new;
v = v_new;
p = p_new;
iterations = iterations + 1;
end
% After the converged solution, we map the staggered variables to
% collocated variables
for i = 1:n_points
for j = 1:n_points
u_final(i,j) = 0.5*(u(i,j) + u(i+1,j));
v_final(i,j) = 0.5*(v(i,j) + v(i,j+1));
p_final(i,j) = 0.25*(p(i,j) + p(i,j+1) + p(i+1,j) + p(i+1,j+1));
end
end
%% Centreline u validation - Comparison with benchmark
figure(11);
plot(u_final(:,(n_points+1)/2),1-y, 'LineWidth', 1)
tab_data = xlsread('plot_u_y_Ghia100.csv', 'A2:B18');
y_ghia = tab_data(:,1);
u_ghia = tab_data(:,2);
figure(11); hold on
plot(u_ghia, y_ghia, 'o', 'LineWidth', 1)
xlabel('u')
ylabel('y')
legend('Numerical', 'Benchmark', 'location', 'southeast')
%% Contour and vector visuals.
x_dom = ((1:n_points)-1).*h;
y_dom = 1-((1:n_points)-1).*h;
[X,Y] = meshgrid(x_dom,y_dom);
figure(21);
contourf(X,Y,u_final, 21, 'LineStyle', 'none')
colorbar
colormap('jet')
xlabel('x')
ylabel('y')
figure(22);
hold on
quiver(X, Y, u_final, v_final, 5, 'k')
%% CASE 10:
case 10
%% Defining the mesh
n_points = 31;
dom_size = 1;
h = dom_size/(n_points - 1);
dt = 0.0001;
alpha = dt/(h*h);
%% Initialising the problem
y(n_points, n_points) = 0;
y(1,:) = 1;
y_new(n_points, n_points) = 0;
y_new(1,:) = 1;
y_transient = 0;
error_mag = 1;
error_req = 1e-6;
iterations = 0;
%Tracking the error magnitude
error_track = 0;
%% Calculations
while error_mag > error_req
for i = 2:(n_points-1)
for j = 2:(n_points-1)
y_new(i,j) = y(i,j) + alpha.*(y(i-1,j) + y(i+1,j) + y(i,j-1) + y(i,j+1) - 4*y(i,j));
y_transient(iterations+1,1:n_points,1:n_points) = y_new;
end
end
iterations = iterations + 1;
% Calculation of error magnitude
error_mag = 0;
for i = 2:(n_points-1)
for j = 2:(n_points-1)
error_mag = error_mag + abs(y(i,j) - y_new(i,j));
error_track(iterations) = error_mag;
end
end
if rem(iterations, 100) == 0
iterations
error_mag
end
% Assigning new to be old
y = y_new;
end
%% Plotting
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,y, 12)
colorbar
% Plot the error with time
figure;
time = dt.*(1:iterations);
plot(time, error_track)
%% Concept of subplots
figure;
subplot(2,1,1)
plot(time, error_track)
subplot(2,1,2)
contourf(X,Y,y, 12)
colorbar
%% Plotting a particular timestep
timestep_selected = 100;
x_dom = ((1:n_points)-1).*h;
y_dom = 1-((1:n_points)-1).*h;
[X,Y] = meshgrid(x_dom,y_dom);
y_timestep = y_transient(timestep_selected, :, :);
y_timestep = reshape(y_timestep, [n_points, n_points]);
contourf(X,Y,y_timestep, 12)
colorbar
title(['Time = ' num2str(timestep_selected*dt) 's'])
%% Animation after every N timesteps
N = 100;
timestep_array = 1:N:iterations;
figure;
for i = 1:length(timestep_array)
timestep_selected = timestep_array(i);
y_timestep = y_transient(timestep_selected, :, :);
y_timestep = reshape(y_timestep, [n_points, n_points]);
contourf(X,Y,y_timestep, 12)
colorbar
title(['Time = ' num2str(timestep_selected*dt) 's'])
pause(0.25);
end
%% CASE 11:
case 11
%% 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
otherwise
disp('input right value');
end