solution for steady flow at Re = 10. in coursera i took mathematics for engineers.for that i need to solve an assignment. For this code i cant get correct stream & vorticity

10 views (last 30 days)
flow_around_cylinder_steady
psi = 101x101
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
omega = 101x101
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Unrecognized function or variable 'plot_Re10'.

Error in solution>flow_around_cylinder_steady (line 49)
plot_Re10(psi);
function [psi, omega] = flow_around_cylinder_steady
Re=10;
%%%%% define the grid %%%%%
n=101; m=101; % number of grid points
N=n-1; M=m-1; % number of grid intervals
h=pi/M; % grid spacing based on theta variable
xi=(0:N)*h; theta=(0:M)*h; % xi and theta variables on the grid
%%%%% Initialize the flow fields %%%%%
psi=zeros(n,m);
omega=zeros(n,m);
psi(n,:)=0 % Free stream boundary condition (psi = 0 at the top edge)
%%%%% Set relax params, tol, extra variables %%%%%
r_psi=1.8; % Set the relaxation parameter here, psi equation
r_omega=0.9; % Set the relaxation parameter here, omega equation
delta=1.e-08; % error tolerance
error=2*delta; % initialize error variable
%%%%% Add any additional variable definitions here %%%%%
...
...
%%%%% Main SOR Loop %%%%%
while (error > delta)
psi_old = psi; omega_old = omega;
for i=2:n-1
for j=2:m-1
psi(i,j)=... % (1 - r_psi) * psi_old(i, j) + ...
r_psi * 0.5 * (psi_old(i+1, j) + psi_old(i-1, j) + ...
psi_old(i, j+1) + psi_old(i, j-1) - h^2 * omega(i, j));
end
end
error_psi=max(abs(psi(:)-psi_old(:)));
omega(1,:)=0 % Boundary condition at theta = 0
for i=2:n-1
for j=2:m-1
omega(i,j)=... % (1 - r_omega) * omega_old(i, j) + ...
r_omega * ( (omega_old(i+1, j) + omega_old(i-1, j) + ...
omega_old(i, j+1) + omega_old(i, j-1) - ...
(4 - (2/h^2)) * omega_old(i, j)) / (4 - (2/h^2)) );
end
end
error_omega=max(abs(omega(:)-omega_old(:)));
error=max(error_psi, error_omega);
end
plot_Re10(psi);
end
Please refer to page 59 (64/69) of the PDF document for the 'plot_Re10' custom plotting function.
  12 Comments
Torsten
Torsten on 8 Aug 2024
And where do you take care of the inscribed cylinder in your code ? I only see a rectangular region in which you compute "psi" and "omega".

Sign in to comment.

Answers (3)

LeoAiE
LeoAiE on 7 Aug 2024
Edited: Sam Chak on 8 Aug 2024
Hi there!
This is not exact solution you looking for but more to help you think about the issue you are facing!
[psi, omega] = flow_around_cylinder_steady
Warning: Contour not rendered for constant ZData
psi = 101x101
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
omega = 101x101
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
function [psi, omega] = flow_around_cylinder_steady
Re = 10;
% Define the grid
n = 101;
m = 101;
N = n - 1;
M = m - 1;
h = 2 * pi / M;
xi = linspace(0, 1, N+1);
theta = linspace(0, 2*pi, M+1);
% Initialize the flow fields
psi = zeros(n, m);
omega = zeros(n, m);
% Boundary conditions for psi (streamfunction)
psi(n, :) = 0; % Top boundary (free stream)
psi(1, :) = 0; % Bottom boundary (cylinder surface)
psi(:, 1) = psi(:, 2); % Left boundary (symmetry)
psi(:, end) = psi(:, end-1); % Right boundary (symmetry)
% Boundary conditions for omega (vorticity)
omega(1, :) = -2 * psi(2, :) / h^2; % No-slip condition at cylinder surface
omega(:, 1) = omega(:, 2); % Left boundary (symmetry)
omega(:, end) = omega(:, end-1); % Right boundary (symmetry)
% Set relax params, tol, extra variables
r_psi = 1.8;
r_omega = 0.9;
delta = 1.e-08;
error = 2 * delta;
% Main SOR Loop
while (error > delta)
psi_old = psi;
omega_old = omega;
% Update psi
for i = 2:n-1
for j = 2:m-1
psi(i, j) = (1 - r_psi) * psi_old(i, j) + ...
r_psi * 0.25 * (psi_old(i+1, j) + psi_old(i-1, j) + ...
psi_old(i, j+1) + psi_old(i, j-1) - h^2 * omega(i, j));
end
end
error_psi = max(abs(psi(:) - psi_old(:)));
% Update omega
for i = 2:n-1
for j = 2:m-1
omega(i, j) = (1 - r_omega) * omega_old(i, j) + ...
r_omega * 0.25 * (omega_old(i+1, j) + omega_old(i-1, j) + ...
omega_old(i, j+1) + omega_old(i, j-1) - ...
(4 - h^2) * omega_old(i, j)) / (4 - h^2);
end
end
error_omega = max(abs(omega(:) - omega_old(:)));
error = max(error_psi, error_omega);
end
plot_Re10(psi);
end
%% Code snippet for plotting 'psi'
function plot_Re10(psi)
% This function plots the streamfunction psi
figure;
contourf(psi', 20); % Transpose psi to get correct orientation
colorbar;
title('Streamfunction \psi for Re = 10');
xlabel('x');
ylabel('y');
end

Sam Chak
Sam Chak on 8 Aug 2024
The code for plot_Re10() function is based on the video lecture of Prof. Jeffrey Chasnov.
[psi, omega] = flow_around_cylinder_steady;
function [psi, omega] = flow_around_cylinder_steady
Re=10;
%%%%% define the grid %%%%%
n=101; m=101; % number of grid points
N=n-1; M=m-1; % number of grid intervals
h=pi/M; % grid spacing based on theta variable
xi=(0:N)*h; theta=(0:M)*h; % xi and theta variables on the grid
%%%%% Initialize the flow fields %%%%%
psi=zeros(n,m);
omega=zeros(n,m);
psi(n,:)=exp(xi(n))*sin(theta(:)); % Write the free stream bc here
%%%%% Set relax params, tol, extra variables %%%%%
r_psi=1.8; % Set the relaxation parameter here, psi equation
r_omega=0.9; % Set the relaxation parameter here, omega equation
delta=1.e-08; % error tolerance
error=2*delta; % initialize error variable
%%%%% Add any additional variable definitions here %%%%%
...
...
%%%%% Main SOR Loop %%%%%
while (error > delta)
psi_old = psi; omega_old = omega;
for i=2:n-1
for j=2:m-1
psi(i,j)=(1-r_psi)*psi(i,j)+(r_psi/4)*(psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1)+h*h*exp(2*xi(i))*omega(i,j));% Write psi equation here
end
end
error_psi=max(abs(psi(:)-psi_old(:)));
omega(1,:)=(psi(3,:) - 8*psi(2,:))/(2*h^2); % Write the boundary condition here
for i=2:n-1
for j=2:m-1
omega(i,j)=(1-r_omega)*omega(i,j)+(r_omega/4)*(omega(i+1,j)+omega(i-1,j)+omega(i,j+1)+omega(i,j-1)+(Re/8)*((psi(i+1,j)-psi(i-1,j))*(omega(i,j+1)-omega(i,j-1))-(psi(i,j+1)-psi(i,j-1))*(omega(i+1,j)-omega(i-1,j)))); % Write omega equation here
end
end
error_omega=max(abs(omega(:)-omega_old(:)));
error=max(error_psi, error_omega);
end
%psi
plot_Re10(psi);
% figure(1)
% contourf(xi,theta,psi.')
% colorbar
% figure(2)
% contourf(xi,theta,omega.')
% colorbar
end
%% Code by Prof. Jeffrey Chasnov
function plot_Re10(psi)
% Reynolds number (design parameter)
Re = 10;
% setup the xi-theta grid
n = size(psi,1);
m = size(psi,2);
N = n - 1;
M = m - 1;
h = pi/M; % grid spacing
xi = (0:N)*h;
theta = (0:M)*h;
[XI, THETA] = meshgrid(xi,theta);
... (Code is very long. Please refer to Prof. Jeffrey Chasnov's video lecture)
end
  6 Comments
Sam Chak
Sam Chak on 9 Aug 2024
This is probably a technical issue of the grader. Perhaps, inform your Professor about issue so that he can liaise with the Coursera technical team.

Sign in to comment.


Cris LaPierre
Cris LaPierre on 9 Aug 2024
Assuming this exercise is coming from Week 2 of Mathematics for Engineers: The Capstone Course, I can't duplicate the error. I'd suggest trying to submit your solution again.

Products


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!