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

Question

SANTHOSH 2024-8-6，9:57

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2143356-solution-for-steady-flow-at-re-10-in-coursera-i-took-mathematics-for-engineers-for-that-i-need-to

回答： Cris LaPierre 2024-8-9，13: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.

https://www.math.hkust.edu.hk/~machas/flow-around-a-cylinder.pdf

12 个评论
显示 10更早的评论隐藏 10更早的评论

SANTHOSH 2024-8-8，10:07

编辑：Torsten 2024-8-8，12:04

在 MATLAB Online 中打开

assignmnet solution.JPG

can anyone pls help me on this. i got the plot, but still i couldnt be able to submit. here is the complete code which i have used. @LeoAiE @Sam Chak @Bhavesh Gyanchandan. i culdnt be able to submit it. it shows "Variable psi has an incorrect value and Variable omega has an incorrect value". i have attached my screen thank you in advance

[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

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

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

Torsten 2024-8-8，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".

SANTHOSH 2024-8-8，15:33

@Torsten. so how to overcome this, could you please help to sort out this

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

LeoAiE 2024-8-7，14:31

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2143356-solution-for-steady-flow-at-re-10-in-coursera-i-took-mathematics-for-engineers-for-that-i-need-to#answer_1495989

编辑：Sam Chak 2024-8-8，8:02

在 MATLAB Online 中打开

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>

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>

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

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

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Sam Chak 2024-8-8，18:33

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2143356-solution-for-steady-flow-at-re-10-in-coursera-i-took-mathematics-for-engineers-for-that-i-need-to#answer_1496754

在 MATLAB Online 中打开

Hi @SANTHOSH

The code for plot_Re10() function is based on the video lecture of Prof. Jeffrey Chasnov.

https://www.youtube.com/watch?v=oTyLgQX7sRU

[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

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

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 个评论
显示 4更早的评论隐藏 4更早的评论

SANTHOSH 2024-8-9，4:29

移动：Sam Chak 2024-8-9，7:11

在 MATLAB Online 中打开

function plot_Re10(psi)
Re=10; 
%%%%% 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);
%%%%% x-y grid %%%%%
nx=640; ny=480/2; % number of pixels in x and half of y
xmin=-1.5; xmax=2.5; ymax=(xmax-xmin)*ny/nx; ymin=-ymax;
x=linspace(xmin,xmax,nx+1); y=linspace(0,ymax,ny+1);
[X,Y]=meshgrid(x,y);
%%%%% construct interpolation points %%%%%
xi_i=0.5*log(x.^2+y.^2);
theta_i=wrapTo2Pi(atan2(Y,X));
%%%%% interpolate %%%%%
psi_xy=interp2(XI,THETA,psi',xi_i,theta_i);
%%%%% set psi to zero inside cylinder %%%%%
psi_xy(xi_i<0)=0;
%%%%% scale contour levels %%%%%
%%%%% negative values have same range as positive values %%%%%
psi_min=min(psi_xy(:));
psi_max=max(psi_xy(:));
psi_xy(psi_xy<0)=psi_xy(psi_xy<0)/abs(psi_min);
psi_xy(psi_xy>0)=psi_xy(psi_xy>0)/abs(psi_max);
%%%%% set colormap for contours %%%%%
levels=linspace(-1,1,1000);
cmap=flipud(jet(length(levels)));
colormap(cmap);
%%%%% plot color contours %%%%%
imagesc(x,y,psi_xy); hold on;
imagesc(x,-y,-psi_xy); % negative values of y
%%%%% plot contour lines %%%%%
v=[-0.9:0.2:-0.1,0,0.0005,0.001,0.002:0.004:0.01,0.02:0.04:0.9];
contour(X,Y,psi_xy,v,'LineColor','k');
contour(X,-Y,-psi_xy,-v,'LineColor','k');
%%%%% draw balck circle for cylinder %%%%%
t=linspace(0,2*pi,1000);
a=cos(t); b=sin(t);
fill(a,b,[0 0 0]);
%%%%% neaten plot %%%%%
set(gca,'YDIR','normal');
axis([xmin xmax ymin ymax]); set(gcf,'color','w'); axis off; axis equal;
text(xmin+0.75*(xmax-xmin),ymin+0.08*(ymax-ymin),...
    ['Re=',num2str(Re,'%3.0f')],'FontSize',22,'color','k');

SANTHOSH 2024-8-9，4:31

移动：Sam Chak 2024-8-9，7:11

在 MATLAB Online 中打开

i followed the thing by prof.jeff, in output it shows not enough input arguments. can anyone pl pls give me code to solve tis question, i am stucked here for past one month, really i dont know ho to solve this, so anyone pls give me correct code for this @Sam Chak @LeoAiE

function plot_Re10(psi)
Re=10; 
%%%%% 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);
%%%%% x-y grid %%%%%
nx=640; ny=480/2; % number of pixels in x and half of y
xmin=-1.5; xmax=2.5; ymax=(xmax-xmin)*ny/nx; ymin=-ymax;
x=linspace(xmin,xmax,nx+1); y=linspace(0,ymax,ny+1);
[X,Y]=meshgrid(x,y);
%%%%% construct interpolation points %%%%%
xi_i=0.5*log(x.^2+y.^2);
theta_i=wrapTo2Pi(atan2(Y,X));
%%%%% interpolate %%%%%
psi_xy=interp2(XI,THETA,psi',xi_i,theta_i);
%%%%% set psi to zero inside cylinder %%%%%
psi_xy(xi_i<0)=0;
%%%%% scale contour levels %%%%%
%%%%% negative values have same range as positive values %%%%%
psi_min=min(psi_xy(:));
psi_max=max(psi_xy(:));
psi_xy(psi_xy<0)=psi_xy(psi_xy<0)/abs(psi_min);
psi_xy(psi_xy>0)=psi_xy(psi_xy>0)/abs(psi_max);
%%%%% set colormap for contours %%%%%
levels=linspace(-1,1,1000);
cmap=flipud(jet(length(levels)));
colormap(cmap);
%%%%% plot color contours %%%%%
imagesc(x,y,psi_xy); hold on;
imagesc(x,-y,-psi_xy); % negative values of y
%%%%% plot contour lines %%%%%
v=[-0.9:0.2:-0.1,0,0.0005,0.001,0.002:0.004:0.01,0.02:0.04:0.9];
contour(X,Y,psi_xy,v,'LineColor','k');
contour(X,-Y,-psi_xy,-v,'LineColor','k');
%%%%% draw balck circle for cylinder %%%%%
t=linspace(0,2*pi,1000);
a=cos(t); b=sin(t);
fill(a,b,[0 0 0]);
%%%%% neaten plot %%%%%
set(gca,'YDIR','normal');
axis([xmin xmax ymin ymax]); set(gcf,'color','w'); axis off; axis equal;
text(xmin+0.75*(xmax-xmin),ymin+0.08*(ymax-ymin),...
    ['Re=',num2str(Re,'%3.0f')],'FontSize',22,'color','k');

SANTHOSH 2024-8-9，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 2024-8-9，7:53

Hi @SANTHOSH

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.

请先登录，再进行评论。

Answer 3

Cris LaPierre 2024-8-9，13:30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2143356-solution-for-steady-flow-at-re-10-in-coursera-i-took-mathematics-for-engineers-for-that-i-need-to#answer_1497184

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.