i have a code for solving the stream function using successive line over relaxation , i should get a graph as metioned below can anyone do correcction tions in this code
13 views (last 30 days)
Show older comments
nx = 41;
ny = 41;
Lx = 1.0;
Ly = 1.0;
A = Lx;
dx = Lx / (nx - 1);
dy = Ly / (ny - 1);
omega = 0.5;
tol = 1e-5;
max_iter = 5000;
% Initialize Psi, zeta, and temperature T
Psi = zeros(nx, ny); % Initial guess for Psi (stream function)
zeta = zeros(nx, ny); % Initial guess for vorticity
T = zeros(nx, ny); % Temperature
% Boundary conditions for initial time tau = 0
T(:, :) = 0; % Initially, T = 0 throughout
zeta(:, :) = 0; % Initially, zeta = 0 throughout
% Boundary conditions for tau > 0
% X = 0 and X = A (left and right boundaries)
Psi(:, 1) = 0; % Psi = 0 at X = 0
Psi(:, end) = 0; % Psi = 0 at X = A
% dPsi/dX = 0 at X = 0 and X = A (already satisfied as Psi = 0 at boundaries)
% dT/dX = 0 at X = 0 and X = A (no-flux boundary)
% Y = 0 (bottom boundary)
Psi(1, :) = 0; % dPsi/dY = 0 at Y = 0
T(1, :) = -1; % T = -1 at Y = 0
% Y = 1 (top boundary)
Psi(end, :) = 0; % dPsi/dY = 0 at Y = 1
T(end, :) = 1; % T = +1 at Y = 1
% Main SLOR iteration loop
for iter = 1:max_iter
Psi_old = Psi; % Save old Psi for comparison
for i = 2:nx-1
for j = 2:ny-1
% Compute the velocity components U and V using centered differences
U = (Psi(i, j+1) - Psi(i, j-1)) / (2 * dy); % U = dPsi/dY
V = -(Psi(i+1, j) - Psi(i-1, j)) / (2 * dx); % V = -dPsi/dX
% Successive Line Over Relaxation formula for Psi
Psi(i, j) = (1 - omega) * Psi(i, j) + ...
omega / 4 * (Psi(i+1, j) + Psi(i-1, j) + ...
Psi(i, j+1) + Psi(i, j-1) + dx^2 * zeta(i, j));
end
end
% Compute the vorticity zeta using centered differences after Psi update
for i = 2:nx-1
for j = 2:ny-1
% zeta = dV/dX - dU/dY
zeta(i, j) = ((Psi(i+1, j) - Psi(i-1, j)) / (2 * dx)) - ...
((Psi(i, j+1) - Psi(i, j-1)) / (2 * dy));
end
end
% Check for convergence using the relative error
error = max(max(abs(Psi - Psi_old)));
if error < tol
fprintf('Converged after %d iterations with error %e\n', iter, error);
break;
end
end
% If maximum iterations are reached without convergence
if iter == max_iter
disp('Maximum iterations reached without convergence.');
end
% Plotting the stream function Psi
figure;
contourf(linspace(0, Lx, nx), linspace(0, Ly, ny), Psi', 20);
colorbar;
title('Stream Function \Psi using SLOR');
xlabel('x');
ylabel('y');
0 Comments
Answers (1)
Malay Agarwal
on 25 Sep 2024
Edited: Malay Agarwal
on 25 Sep 2024
The problem is with the initial guess for Ψ. Since you initialize it as zeros, the following equation:
(1 - omega) * Psi(i, j) + ...
omega / 4 * (Psi(i+1, j) + Psi(i-1, j) + ...
Psi(i, j+1) + Psi(i, j-1) + dx^2 * zeta(i, j))
Always yields zero. Thus, Psi is never updated and the SOR converges after a single iteration with all values in Psi set to zero. This leads to the contour plot not being as expected. I am not sure how Psi should be initialized but if I just do random initialization (refer to the attached file), I do get something sensible:
sor
Note that I have also changed the number of levels to 7 since that's how many I counted in your expected plot.
Hope this helps!
4 Comments
Malay Agarwal
on 25 Sep 2024
Please check at the top of the answer, just above the "Hi". There's a file called "sor.m".
See Also
Categories
Find more on Computational Fluid Dynamics (CFD) in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!