Filter löschen
Filter löschen

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

113 Ansichten (letzte 30 Tage)
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 Kommentare
Torsten
Torsten am 8 Aug. 2024 um 15:15
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".

Melden Sie sich an, um zu kommentieren.

Antworten (3)

LeoAiE
LeoAiE am 7 Aug. 2024 um 14:31
Bearbeitet: Sam Chak am 8 Aug. 2024 um 8:02
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 am 8 Aug. 2024 um 18:33
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 Kommentare
SANTHOSH
SANTHOSH am 9 Aug. 2024 um 7:31
hi @Sam Chak thank you so much for your time and consideration.
but how can i submit it in matrix form
Sam Chak
Sam Chak am 9 Aug. 2024 um 7:53
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.

Melden Sie sich an, um zu kommentieren.


Cris LaPierre
Cris LaPierre am 9 Aug. 2024 um 13:30
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.

Kategorien

Mehr zu Contour Plots finden Sie in Help Center und File Exchange

Produkte


Version

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by