Mathematics for Engineers: The Capstone Course, Assignment 2

omega = flow_around_cylinder_unsteady

Unrecognized function or variable 'plot_Re60'.

Error in solution>flow_around_cylinder_unsteady (line 65)
plot_Re60(omega,t_end);

function omega = flow_around_cylinder_unsteady
Re=60;
%%%%% define the grid %%%%%
n=101; m=202; % number of grid points
N=n-1; M=m-2; % number of grid intervals: 2 ghost points, theta=-h,2*pi
h=2*pi/M; % grid spacing based on theta variable
xi=(0:N)*h; theta=(-1:M)*h; % xi and theta variables on the grid
%%%%% time-stepping parameters %%%%%
t_start=0; t_end=0.5; %vortex street starts at around t=1000
tspan=[t_start t_end]; 
%%%%% Initialize vorticity field %%%%%
omega=zeros(n,m);
%%%%% Construct the matrix A for psi equation %%%%%
boundary_index_bottom = 1:n;
boundary_index_left = 1:n:1+(m-1)*n;
boundary_index_top = 1+(m-1)*n:m*n;
boundary_index_right = n:n:m*n;
diagonals1 = [4*ones(n*m,1), -ones(n*m,4)];
A=spdiags(diagonals1,[0 -1 1 -n n], n*m, n*m);
I=speye(m*n);
A(boundary_index_left,:)=I(boundary_index_left,:);
A(boundary_index_right,:)=I(boundary_index_right,:);
diagonals2 = [ones(n*m,1), -ones(n*m,2)];
Aux = spdiags(diagonals2,[0 n -n],n*m,n*m);
A(boundary_index_top,:)=Aux(boundary_index_top,:);
A(boundary_index_bottom,:)=Aux(boundary_index_bottom,:);

%%%%% Find the LU decomposition %%%%%
[L,U]=lu(A); clear A;
%%%%% Compute any time-independent constants %%%%%

%%%%% advance solution using ode23 %%%%%
options=odeset('RelTol', 1.e-03);
omega=omega(2:n-1,2:m-1); % strip boundary values for ode23
omega=omega(:); % make a column vector
[t,omega]=ode23...
    (@(t,omega)omega_eq(omega,L,U, theta, xi, h, Re, n, m),...
    tspan, omega, options);
%%%%% expand omega to include boundaries %%%%%
temp=zeros(n,m);
temp(2:n-1,2:m-1)=reshape(omega(end,:),n-2,m-2);
omega=temp; clear temp;
%%%%% compute stream function (needed for omega boundary values) %%%%%
omega_tilde = zeros(n, m);
omega_tilde(1, :) = 0;
omega_tilde(n, :) = exp(xi(n)) * sin(theta);
omega_tilde(:, 1) = 0;
omega_tilde(:, m) = 0;
for i = 2:n-1
    for j = 2:m-1
        omega_tilde(i, j) = h^2 * exp(2 * xi(i)) * omega(i, j);
    end
end
omega_tilde = omega_tilde(:);
psi = reshape(U \ (L \ omega_tilde), n, m);

%%%%% set omega boundary conditions %%%%%
omega(1, :) = (psi(3, :) - 8 * psi(2, :)) * 1 / (2 * h^2);
omega(n, :) = 0;
omega(:, 1) = omega(:, m-1);
omega(:, m) = omega(:, 2);
%%%%% plot scalar vorticity field %%%%%
plot_Re60(omega,t_end);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d_omega_dt=omega_eq(omega,L,U,theta, xi, h, Re, n, m)
%%%%% expand omega to include boundary points %%%%%
temp=zeros(n,m);
index1=2:n-1; index2=2:m-1;
temp(index1,index2)=reshape(omega,n-2,m-2);
omega=temp; clear temp;
%%%%% compute stream function %%%%%
omega_tilde = zeros(n, m);
omega_tilde(1, :) = 0;
omega_tilde(n, :) = exp(xi(n)) * sin(theta);
omega_tilde(:, 1) = 0;
omega_tilde(:, m) = 0;
for i = 2:n-1
    for j = 2:m-1
        omega_tilde(i, j) = h^2 * exp(2 * xi(i)) * omega(i, j);
    end
end
omega_tilde = omega_tilde(:);
psi = reshape(U \ (L \ omega_tilde), n, m);
%%%%% compute derivatives of omega %%%%%
omega(1, :) = (psi(3, :) - 8 * psi(2, :)) * 1 / (2 * h^2);
omega(n, :) = 0;
omega(:, 1) = omega(:, m-1);
omega(:, m) = omega(:, 2);
d_omega_dt = zeros(n-2, m-2);
for i = 2:n-1
    for j = 2:m-1
        d_omega_dt(i-1, j-1) = ...
            2 * exp(-2 * xi(i)) * 1 / (h^2 * Re) * ...
            (omega(i+1, j) + omega(i-1, j) + omega(i, j+1) + omega(i, j-1) - 4 * omega(i, j)) ...
            + exp(-2 * xi(i)) / (4 * h^2) * ...
            ((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)));
    end
end
d_omega_dt = d_omega_dt(:);
end

[Ayush]

Can anyone find the error in my submission

[Ayush]

the grey line cannot be changed and plot function is defined somewhere else we cannot change that too

[Ayush]

[Ayush]

[Walter Roberson]

Is plot_Re60 the same code as is posted at https://www.mathworks.com/matlabcentral/answers/1815940-can-anyone-please-help-me-with-this#comment_2472053 ?

[Ayush]

@Walter Roberson could you please provide me the solution of this code

[Ayush]

@Walter Roberson yes it is the same

[Walter Roberson]

You have not indicated what problem you are encountering.

Your code runs and produces something so it is not immediately clear that anything is wrong at all.

[Ayush]

yes it does produces the result but the problem is while submitting the evaluator is showing inncorrect values.

[Torsten]

Is the problem similar to the one reported here:

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

?

[Ayush]

@Torsten no

[Walter Roberson]

Well, there are several different possibilities:

• your equations might be wrong
• your equations might be right, but due to round-off effects the results are being determined to be incorrect
• your equations might be right and the values might be right to within configured round-off, but there might be a glitch in approving the results (for example, the internal test might have been configured as ~= instead of ==)