can anyone please help me with this ?

im trying to compute the steady solutions for the stream function and scalar vorticity for a 2d flow around an infinite cylinder at RE = 10 Here is the code:

[psi, omega] = flow_around_cylinder_steady

psi = 101×101

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0010 0.0020 0.0030 0.0039 0.0049 0.0059 0.0069 0.0078 0.0088 0.0097 0.0106 0.0116 0.0125 0.0134 0.0143 0.0151 0.0160 0.0168 0.0177 0.0185 0.0193 0.0200 0.0208 0.0215 0.0222 0.0229 0.0236 0.0242 0.0248 0 0.0020 0.0039 0.0059 0.0079 0.0098 0.0118 0.0137 0.0156 0.0175 0.0194 0.0213 0.0231 0.0249 0.0267 0.0285 0.0302 0.0320 0.0336 0.0353 0.0369 0.0385 0.0400 0.0415 0.0430 0.0444 0.0458 0.0471 0.0484 0.0496 0 0.0030 0.0059 0.0089 0.0118 0.0147 0.0176 0.0205 0.0234 0.0263 0.0291 0.0319 0.0346 0.0374 0.0401 0.0427 0.0453 0.0479 0.0504 0.0529 0.0553 0.0577 0.0600 0.0622 0.0644 0.0665 0.0686 0.0706 0.0725 0.0744 0 0.0039 0.0079 0.0118 0.0157 0.0196 0.0235 0.0273 0.0312 0.0350 0.0387 0.0425 0.0461 0.0498 0.0534 0.0569 0.0604 0.0638 0.0672 0.0704 0.0737 0.0768 0.0799 0.0829 0.0858 0.0886 0.0914 0.0940 0.0966 0.0990 0 0.0049 0.0098 0.0147 0.0196 0.0245 0.0293 0.0341 0.0389 0.0436 0.0483 0.0530 0.0576 0.0621 0.0666 0.0710 0.0754 0.0796 0.0838 0.0879 0.0919 0.0959 0.0997 0.1035 0.1071 0.1106 0.1140 0.1173 0.1205 0.1236 0 0.0059 0.0118 0.0176 0.0235 0.0293 0.0351 0.0409 0.0466 0.0523 0.0579 0.0635 0.0690 0.0744 0.0798 0.0851 0.0903 0.0954 0.1004 0.1053 0.1101 0.1148 0.1194 0.1239 0.1283 0.1325 0.1366 0.1406 0.1444 0.1481 0 0.0069 0.0137 0.0205 0.0273 0.0341 0.0409 0.0476 0.0543 0.0609 0.0674 0.0739 0.0803 0.0866 0.0929 0.0990 0.1051 0.1110 0.1169 0.1226 0.1282 0.1337 0.1390 0.1443 0.1493 0.1543 0.1590 0.1636 0.1681 0.1724 0 0.0078 0.0156 0.0234 0.0312 0.0389 0.0466 0.0543 0.0618 0.0694 0.0768 0.0842 0.0915 0.0988 0.1059 0.1129 0.1198 0.1266 0.1333 0.1398 0.1462 0.1524 0.1585 0.1645 0.1702 0.1759 0.1813 0.1865 0.1916 0.1965 0 0.0088 0.0175 0.0263 0.0350 0.0436 0.0523 0.0609 0.0694 0.0778 0.0862 0.0945 0.1027 0.1108 0.1188 0.1267 0.1344 0.1420 0.1495 0.1568 0.1640 0.1710 0.1778 0.1845 0.1910 0.1973 0.2034 0.2093 0.2150 0.2204

omega = 101×101

-2.9995 -2.9995 -5.9961 -8.9867 -11.9684 -14.9384 -17.8936 -20.8311 -23.7481 -26.6417 -29.5089 -32.3471 -35.1533 -37.9248 -40.6589 -43.3529 -46.0041 -48.6099 -51.1677 -53.6750 -56.1294 -58.5283 -60.8695 -63.1506 -65.3694 -67.5237 -69.6114 -71.6303 -73.5786 -75.4542 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

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(:));
%%%%% Set relax params, tol, extra variables %%%%%
r_psi=1.8;
r_omega=0.9;
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)=sin(xi(i)) * sin(theta(j));
        end
    end
    error_psi=max(abs(psi(:)-psi_old(:)));
    omega(1,:)= (psi(3,j) - 8*psi(2,j)) * 1/(2*h^2);
    for i=2:n-1
        for j=2:m-1
            omega(1,j)= (psi(3,j) - 8*psi(2,j)) * 1/(2*h^2);
        end
    end
    error_omega=max(abs(omega(:)-omega_old(:)));
    error=max(error_psi, error_omega);
end
%plot_Re10(psi);
end

The code to call the function is: [psi, omega] = flow_around_cylinder_steady;

The gray areas cannot be changed

but i keep getting a message saying that my psi and omega variables are incorrect.

can anyone please help me pinpoint the exact area of my problem and let me know how to fix it ?

[Ghazwan]

can you please format your question in a proper way?

https://www.mathworks.com/matlabcentral/answers/7885-tutorial-how-to-format-your-question

[Torsten]

@Mohamed

As you can see above, no such message pops up. Where is the function plot_Re10 ?

[Mohamed]

the answer im getting looks like this:

but the actual answer looks like this:

the gray areas in the question cannot be changed

[SANTHOSH]

but after run the function i couldnt be able to submit. it shows that " Variable psi has an incorrect value and Variable omega has an incorrect value"

[SANTHOSH]

@Mohamed Mohamed @Torsten 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. 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

Answer 1

[Aisha Singhania]

There are some things you need to change, try some of these codes. Although I am not sure myself but hopefully this does help you. These were some of the errors I had noticed in your code.

%%%%% 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)+fac(i)*omega(i,j));
        end
    end
    error_psi=max(abs(psi(:)-psi_old(:)));
      omega(1,:)=(1/(2*h^2))*(psi(3,:)-8*psi(2,:)); %Second order bc
    for i=2:n-1
        for j=2:m-1
            f(i,j)=...
                (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));
            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)*f(i,j));
        end
    end
    error_omega=max(abs(omega(:)-omega_old(:)));
    error=max(error_psi, error_omega);
end
plot_Re10(psi);