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)
Ältere Kommentare anzeigen
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,:)=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
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".
Antworten (3)
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
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
0 Kommentare
Sam Chak
am 8 Aug. 2024 um 18:33
Hi @SANTHOSH
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
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.
0 Kommentare
Siehe auch
Kategorien
Mehr zu Contour Plots finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!