# 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

71 views (last 30 days)

Show older comments

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 Comments

Torsten
on 8 Aug 2024 at 15:15

### Answers (3)

LeoAiE
on 7 Aug 2024

Edited: Sam Chak
on 8 Aug 2024 at 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 Comments

Sam Chak
on 8 Aug 2024 at 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 Comments

Cris LaPierre
on 9 Aug 2024 at 13:30

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!