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 个评论
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