clear; clc;
% Parameters
n_x = 30;
n_y = 30;
dx = 1.0 / n_x;
dy = 1.0 / n_y;
Re = 100;
% Function to compute momentum link coefficients
function [A_p, A_e, A_w, A_n, A_s, source_x, source_y] =
momentum_link_coefficients(u_star, u_face, v_face, p, source_x, source_y, A_p, A_e,
A_w, A_n, A_s, n_x, n_y, dx, dy, Re)
D_e = dy / (dx * Re);
D_w = dy / (dx * Re);
D_n = dx / (dy * Re);
D_s = dx / (dy * Re);
for i = 2:n_y
for j = 2:n_x
F_e = dy * u_face(i, j);
F_w = dy * u_face(i, j-1);
F_n = dx * v_face(i-1, j);
F_s = dx * v_face(i, j);
A_e(i, j) = D_e + max(0.0, -F_e);
A_w(i, j) = D_w + max(0.0, F_w);
A_n(i, j) = D_n + max(0.0, -F_n);
A_s(i, j) = D_s + max(0.0, F_s);
A_p(i, j) = A_w(i, j) + A_e(i, j) + A_n(i, j) + A_s(i, j) + (F_e - F_w) +
(F_n - F_s);
source_x(i, j) = 0.5 * (p(i, j-1) - p(i, j+1)) * dx;
source_y(i, j) = 0.5 * (p(i+1, j) - p(i-1, j)) * dy;
end
end
% Boundary conditions handling
% Similar boundary condition handling for walls and corners as in Python code
% ...
end
% Function to solve the equations iteratively
function [u, norm] = solve(u, u_star, A_p, A_e, A_w, A_n, A_s, source_x, alpha,
epsilon, max_inner_iteration, l2_norm, n_x, n_y)
for n = 1:max_inner_iteration
l2_norm = 0;
for i = 1:n_y+1
for j = 1:n_x+1
u(i, j) = alpha * (A_e(i, j) * u(i, j+1) + A_w(i, j) * u(i, j-1) +
A_n(i, j) * u(i-1, j) + A_s(i, j) * u(i+1, j) + source_x(i, j)) / A_p(i, j) + (1 -
alpha) * u_star(i, j);
l2_norm = l2_norm + (u(i, j) - alpha * (A_e(i, j) * u(i, j+1) +
A_w(i, j) * u(i, j-1) + A_n(i, j) * u(i-1, j) + A_s(i, j) * u(i+1, j) + source_x(i,
j)) / A_p(i, j) - (1 - alpha) * u_star(i, j))^2;
end
end
if n == 1
norm = sqrt(l2_norm);
end
l2_norm = sqrt(l2_norm);
if l2_norm < epsilon
break;
end
end
end
% Function to calculate face velocities
function [u_face, v_face] = face_velocity(u, v, u_face, v_face, p, A_p, alpha_uv,
n_x, n_y, dy)
for i = 1:n_y+1
for j = 1:n_x
u_face(i, j) = 0.5 * (u(i, j) + u(i, j+1)) + 0.25 * alpha_uv * (p(i, j+1)
- p(i, j-1)) * dy / A_p(i, j) + 0.25 * alpha_uv * (p(i, j+2) - p(i, j)) * dy / A_p(i,
j+1) - 0.5 * alpha_uv * (1 / A_p(i, j) + 1 / A_p(i, j+1)) * (p(i, j+1) - p(i, j)) *
dy;
end
end
for i = 2:n_y+1
for j = 1:n_x+1
v_face(i-1, j) = 0.5 * (v(i, j) + v(i-1, j)) + 0.25 * alpha_uv * (p(i-1,
j) - p(i+1, j)) * dy / A_p(i, j) + 0.25 * alpha_uv * (p(i-2, j) - p(i, j)) * dy /
A_p(i-1, j) - 0.5 * alpha_uv * (1 / A_p(i, j) + 1 / A_p(i-1, j)) * (p(i-1, j) - p(i,
j)) * dy;
end
end
end
% Function to correct the pressure
function p_star = correct_pressure(p_star, p, p_prime, alpha_p, n_x, n_y)
p_star = p + alpha_p * p_prime;
% Boundary conditions
p_star(1, 2:n_x+1) = p_star(2, 2:n_x+1);
p_star(2:n_y+1, 1) = p_star(2:n_y+1, 2);
p_star(2:n_y+1, n_x+2) = p_star(2:n_y+1, n_x+1);
p_star(n_y+2, 2:n_x+1) = p_star(n_y+1, 2:n_x+1);
p_star(1, 1) = (p_star(2, 2) + p_star(1, 2) + p_star(2, 1)) / 3;
p_star(1, n_x+2) = (p_star(1, n_x+1) + p_star(2, n_x+1) + p_star(2, n_x+2)) / 3;
p_star(n_y+2, 1) = (p_star(n_y+1, 1) + p_star(n_y+1, 2) + p_star(n_y+2, 2)) / 3;
p_star(n_y+2, n_x+2) = (p_star(n_y+1, n_x+2) + p_star(n_y+2, n_x+1) +
p_star(n_y+1, n_x+1)) / 3;
end
% Function to correct cell center velocities
function [u_star, v_star] = correct_cell_center_velocity(u, v, u_star, v_star,
p_prime, A_p, alpha_uv, n_x, n_y, dy, dx)
for i = 1:n_y+1
for j = 2:n_x
u_star(i, j) = u(i, j) + 0.5 * alpha_uv * (p_prime(i, j-1) - p_prime(i,
j+1)) * dy / A_p(i, j);
end
end
for i = 2:n_y
for j = 1:n_x+1
v_star(i, j) = v(i, j) + 0.5 * alpha_uv * (p_prime(i+1, j) - p_prime(i-1,
j)) * dx / A_p(i, j);
end
end
end
% Function to correct face velocities
function [u_face, v_face] = correct_face_velocity(u_face, v_face, p_prime, A_p,
alpha_uv, n_x, n_y, dy, dx)
for i = 1:n_y+1
for j = 1:n_x
u_face(i, j) = u_face(i, j) + 0.5 * alpha_uv * (1 / A_p(i, j) + 1 /
A_p(i, j+1)) * (p_prime(i, j) - p_prime(i, j+1)) * dy;
end
end
for i = 2:n_y+1
for j = 1:n_x+1
v_face(i-1, j) = v_face(i-1, j) + 0.5 * alpha_uv * (1 / A_p(i, j) + 1 /
A_p(i-1, j)) * (p_prime(i, j) - p_prime(i-1, j)) * dx;
end
end
end
% Function to visualize results
function post_processing(u_star, v_star, p_star, X, Y, x, y)
figure(1);
contourf(X, Y, flipud(u_star), 50, 'LineColor', 'none');
colorbar;
title('U contours');
figure(2);
contourf(X, Y, flipud(v_star), 50, 'LineColor', 'none');
colorbar;
title('V contours');
figure(3);
contourf(X, Y, flipud(p_star), 50, 'LineColor', 'none');
colorbar;
title('P contours');
figure(4);
plot(1-y, u_star(:, round(n_x / 2)));
xlabel('y');
ylabel('u');
title('U centerline velocity');
figure(5);
plot(x, v_star(round(n_y / 2), :));
xlabel('x');
ylabel('v');
title('V centerline velocity');
end
% Main
u = zeros(n_y+2, n_x+2);
u_star = zeros(n_y+2, n_x+2);
v = zeros(n_y+2, n_x+2);
v_star = zeros(n_y+2, n_x+2);
p_star = zeros(n_y+2, n_x+2);
p = zeros(n_y+2, n_x+2);
p_prime = zeros(n_y+2, n_x+2);
A_p = ones(n_y+2, n_x+2);
A_e = ones(n_y+2, n_x+2);
A_w = ones(n_y+2, n_x+2);
A_s = ones(n_y+2, n_x+2);
A_n = ones(n_y+2, n_x+2);
Ap_p = ones(n_y+2, n_x+2);
Ap_e = ones(n_y+2, n_x+2);
Ap_w = ones(n_y+2, n_x+2);
Ap_s = ones(n_y+2, n_x+2);
Ap_n = ones(n_y+2, n_x+2);
source_x = zeros(n_y+2, n_x+2);
source_y = zeros(n_y+2, n_x+2);
source_p = zeros(n_y+2, n_x+2);
u_face = zeros(n_y+2, n_x+1);
v_face = zeros(n_y+1, n_x+2);
x = [0, linspace(dx/2, 1-dx/2, n_x), 1];
y = [0, linspace(dy/2, 1-dy/2, n_y), 1];
[X, Y] = meshgrid(x, y);
u(1, 2:n_x+1) = 1;
u_star(1, 2:n_x+1) = 1;
u_face(1, 1:n_x) = 1;
l2_norm_x = 0;
alpha_uv = 0.7;
epsilon_uv = 1e-3;
max_inner_iteration_uv = 50;
l2_norm_y = 0;
l2_norm_p = 0;
max_inner_iteration_p = 200;
dummy_alpha_p = 1;
epsilon_p = 1e-4;
alpha_p = 0.2;
max_outer_iteration = 200;
for n = 1:max_outer_iteration
[A_p, A_e, A_w, A_n, A_s, source_x, source_y] =
momentum_link_coefficients(u_star, u_face, v_face, p, source_x, source_y, A_p, A_e,
A_w, A_n, A_s, n_x, n_y, dx, dy, Re);
[u, l2_norm_x] = solve(u, u_star, A_p, A_e, A_w, A_n, A_s, source_x, alpha_uv,
epsilon_uv, max_inner_iteration_uv, l2_norm_x, n_x, n_y);
[v, l2_norm_y] = solve(v, v_star, A_p, A_e, A_w, A_n, A_s, source_y, alpha_uv,
epsilon_uv, max_inner_iteration_uv, l2_norm_y, n_x, n_y);
[u_face, v_face] = face_velocity(u, v, u_face, v_face, p, A_p, alpha_uv, n_x,
n_y, dy);
[Ap_p, Ap_e, Ap_w, Ap_n, Ap_s, source_p] =
pressure_correction_link_coefficients(u, u_face, v_face, Ap_p, Ap_e, Ap_w, Ap_n,
Ap_s, source_p, A_p, A_e, A_w, A_n, A_s, alpha_uv, n_x, n_y, dx, dy);
[p_prime, l2_norm_p] = solve(p_prime, p_prime, Ap_p, Ap_e, Ap_w, Ap_n, Ap_s,
source_p, dummy_alpha_p, epsilon_p, max_inner_iteration_p, l2_norm_p, n_x, n_y);
p_star = correct_pressure(p_star, p, p_prime, alpha_p, n_x, n_y);
[u_star, v_star] = correct_cell_center_velocity(u, v, u_star, v_star, p_prime,
A_p, alpha_uv, n_x, n_y, dy, dx);
[u_face, v_face] = correct_face_velocity(u_face, v_face, p_prime, A_p, alpha_uv,
n_x, n_y, dy, dx);
p = p_star;
disp([l2_norm_x, l2_norm_y, l2_norm_p]);
if l2_norm_x < 1e-4 && l2_norm_y < 1e-4 && l2_norm_p < 1e-4
disp('Converged!');
break;
end
end
post_processing(u_star, v_star, p_star, X, Y, x, y);
