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

조회 수: 29 (최근 30일)
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
Torsten
Torsten 2024년 8월 8일
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".

댓글을 달려면 로그인하십시오.

답변 (3개)

LeoAiE
LeoAiE 2024년 8월 7일
편집: Sam Chak 2024년 8월 8일
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 2024년 8월 8일
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
SANTHOSH
SANTHOSH 2024년 8월 9일
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 2024년 8월 9일
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.

댓글을 달려면 로그인하십시오.


Cris LaPierre
Cris LaPierre 2024년 8월 9일
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.

카테고리

Help CenterFile Exchange에서 Legend에 대해 자세히 알아보기

제품


릴리스

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by