Contour plot_fluid mechanics

Hello, I'm trying to make a 2D contour plot of an advanced fluid mechanical equation. The results I get, are not logical at all and I cannot locate the mistake. The code is:

function [U,L,out15]=contour_plot

Ca=1; Cd=1; Cl=1; z=0.001; m=10; S=0.2;

u=linspace(1, 14);

l=linspace(0.5, 3.16); [U,L]= meshgrid(u,l);

out15= zeros(length(U), length(L));

for i= 1:length(U)

 for j= 1:length(L)

      fprintf('Calculating for %d\n', i, j);

[t,y]=ode23s(@(t,y) motion5(t,y, i, j), [0,100],[0,0]);

equation_of_power(t,y,z,m, i, j);

out15(i,j)= equation_of_power(t,y,z,m, i, j);

 end

end

%figure(1), plot(t,y(:,1)), pause(1)

     function power12= equation_of_power(t,y,z,m, i,j)
 Y = y(:,2);

 equation_of_power=4*pi*z*m*(1./U(i)).^3.*(L(j)^2)*Y.^2;

power12 = trapz(t,equation_of_power);

end

 function dydt=motion5(t,y, i, j)
    k1=m*(L(j)./U(i)).^2+Ca*(L(j)./U(i)).^2*(((((L(j)./U(i)).*y(2))^2+2*(L(j)./U(i)).*y(2)*sin(y(1))+(sin(y(1)))^2)*(1/(1+((L(j)./U(i)).*y(2))^2+2*(L(j)./U(i)).*y(2)*sin(y(1))))));
dydt=[y(2);(1/k1)*((((-4*pi*m*z)*(L(j)./U(i)).^2)-....
    (L(j)./U(i))*Ca*(((L(j)./U(i)).*y(2))^2*cos(y(1))+(L(j)./U(i)).*y(2)*sin(y(1))*cos(y(1)))*(1/(1+((L(j)./U(i)).*y(2))^2+2*(L(j)./U(i)).*y(2)*sin(y(1))))-...
    (2/pi)*(L(j)^2)*(1./U(i)).*Cd*sqrt(1+((L(j)./U(i)).*y(2))^2+2*(L(j)./U(i)).*y(2)*sin(y(1))))*y(2)-((4*(pi)^2)*m*(L(j)./U(i)).^2).*y(1)+...
    ((2/pi)*Cl*L(j)*cos(y(1))*sin(2*pi.*U(i).*S*t)*(1/sqrt((1+((L(j)./U(i)).*y(2))^2+2*(L(j)./U(i)).*y(2)*sin(y(1))))))-...
    ((2/pi)*L(j)*Cd*sin(y(1))*sqrt(1+((L(j)./U(i)).*y(2))^2+2*(L(j)./U(i)).*y(2)*sin(y(1)))))];
end

figure

contourf(U,L,out15);

end

[Tim Berk]

I tried running your code, but gave up after a couple of minutes. You are solving your differential equation 10000 times which is a bit too much for my computer/time.

Are you sure looping is the best way? Can't you linearize the process (i.e. do matrix calculations instead of looping through each individual point of the matrix)?

It would be helpful if you can explain what you expected to see vs what you see.

[Iro Malefaki]

Hello Tim. Thank you for answering. The result I get after running the code, is the attached image. The problem is that output "out15" has very small prices, they should be around 0.15-0.45.

Answer 1

[Tim Berk]

The problem is that you are creating two matrices (U and L) where U has constant columns, while L has constant rows. You loop over i=1:length(U) (which is 1:100) and over j=1:length(L) (also 1:100).

You then call values of L as L(j), which is fine, L(1:100) are the first 100 entries of L, which is the first column i.e. the values in l.

But you also call value of U as U(i), which returns the first 100 values of U, which is the first column, which is the same value 1 each time.

• I'm not sure why you need the meshgrid, but I think you could just loop over i=1:length(u) and j=1:length(l), calling u(i) and l(j) in your ODE.
• Alternatively, call the values as L(1,j) and U(i,1).
• Or loop over the entire meshgrid ( and be prepared to wait a couple of days!!) i.e.i = 1:length(U(:)) and j = 1:length(L(:))