Contour plot_fluid mechanics

4 次查看（过去 30 天）
Iro Malefaki
Iro Malefaki 2017-9-19
评论： 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
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
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
Iro Malefaki 2017-9-19
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.

请先登录，再进行评论。

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

回答（1 个）

Tim Berk
Tim Berk 2017-9-19
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(:))

类别

在 Help Center 和 File Exchange 中查找有关 Fluid Dynamics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by