Contour plot_fluid mechanics

Question

Iro Malefaki 2017-9-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/357375-contour-plot_fluid-mechanics

评论： Iro Malefaki 2017-9-19

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

%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

2 个评论
显示无隐藏无

Tim Berk 2017-9-19

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 2017-9-19

isxusSunathseiUkaiL1.jpg

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 2017-9-19

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/357375-contour-plot_fluid-mechanics#answer_282213

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(:))