Plotting a function given a range of inputs

조회 수: 107 (최근 30일)
Jesse Green
Jesse Green 2018년 3월 6일
댓글: Vatsal Dhamelia 2020년 5월 11일
I need to plot a line of the Moody chart given a relative roughness and range of Reynolds numbers but I have no clue how to write a plot code. For simplicity, lets say I want to plot y=8x+10 with a range of 2<x<20. How would I go about doing this? Thank you.

답변 (2개)

Elias Gule
Elias Gule 2018년 3월 6일
First define your vector x as:
dx = 0.01;
x = 2:dx:20; % where dx in the increment from one x value to the other.
OR
N = 40; % number of samples
x = linspace(2,20,N);
then define y as:
y = 8*x + 10;
Then to plot:
plot(x,y);
For more help please look at the VERY USEFUL Matlab documentation. In the command window just type:
doc plot
  댓글 수: 2
Vatsal Dhamelia
Vatsal Dhamelia 2020년 5월 11일
clear;
clc;
pi = acos(-1);
Um = 1.1; % Fluid Velicity amplitude
D = 0.11; % cylinder diamter
L = 1; % cylinder length
m = 50; % cylinder mass
rho = 1024; % fluid density
K = 200; % spring stiffness
c = 100; % damping constant
CA = 1; % added mass coefficient
CD = 1.8; % drag coefficient
KC = 2:2:20;%KC number
T = KC*D/Um ; % period of oscillatory flow ;
omega = 2*pi/T; % angular frequency of the flow
md = rho*pi*(D*D)/4*L; % added mass coefficient
Ap = D*L; % projected area
dt = T/40; % time step
ndt = T/dt*5; % %total step to be calculated
X(1:ndt+1) = 0; % save displacement X from step 0 to step ndt
V(1:ndt+1) = 0;
Pa(1:ndt+1) = 0;
time(1:ndt+1) = (0:ndt)*dt;
for n=1:ndt
% to calculate kx1 and kv1
ta = time(n);
aw = (omega *Um*cos(omega *ta))-((2*omega*Um/3)*cos(2*omega*ta)); % acceleration of the fluid
Vw = (Um*sin(omega*ta))-((Um/3)*sin(2*omega*ta)); % velocity of the fluid
Va = V(n);
Xa = X(n);
t1 = CA*md*aw;
t2 = 0.5*rho*CD*Ap*abs(Vw-Va)*(Vw-Va);
t3 = -K*Xa-c*Va;
kx1 = V(n);
kv1 = (t1+t2+t3)/(m+CA*md);
% to calculate kx2 and kv2
ta = time(n)+dt/2;
aw = (omega *Um*cos(omega *ta))-((2*omega *Um/3)*cos(2*omega*ta)); % acceleration of the fluid
Vw = Um*sin(omega*ta)-(Um/3)*sin(2*omega*ta); % velocity of the fluid
Va = V(n)+0.5*dt*kv1;
Xa = X(n)+0.5*dt*kx1;
t1 = CA*md*aw;
t2 = 0.5*rho*CD*Ap*abs(Vw-Va)*(Vw-Va);
t3 = -K*Xa-c*Va;
kx2 = V(n)+0.5*dt*kv1;
kv2 = (t1+t2+t3)/(m+CA*md);
% to calculate kx3 and kv3
ta = time(n)+dt/2;
aw = (omega *Um*cos(omega *ta))-((2*omega *Um/3)*cos(2*omega*ta)); % acceleration of the fluid
Vw = Um*sin(omega*ta)-(Um/3)*sin(2*omega*ta); % velocity of the fluid
Va = V(n)+0.5*dt*kv2;
Xa = X(n)+0.5*dt*kx2;
t1 = CA*md*aw;
t2 = 0.5*rho*CD*Ap*abs(Vw-Va)*(Vw-Va);
t3 = -K*Xa-c*Va;
kx3 = V(n)+0.5*dt*kv2;
kv3 = (t1+t2+t3)/(m+CA*md);
% to calculate kx4 and kv4
ta = time(n)+dt;
aw = (omega *Um*cos(omega *ta))-((2*omega *Um/3)*cos(2*omega*ta)); % acceleration of the fluid
Vw = Um*sin(omega*ta)-(Um/3)*sin(2*omega*ta); % velocity of the fluid
Va = V(n)+0.5*dt*kv3;
Xa = X(n)+0.5*dt*kx3;
t1 = CA*md*aw;
t2 = 0.5*rho*CD*Ap*abs(Vw-Va)*(Vw-Va);
t3 = -K*Xa-c*Va;
kx4 = V(n)+dt*kv3;
kv4 = (t1+t2+t3)/(m+CA*md);
% next step values
X(n+1) = X(n)+(dt/6)*(kx1+2*kx2+2*kx3+kx4);
V(n+1) = V(n)+(dt/6)*(kv1+2*kv2+2*kv3+kv4);
Power = (CA*md*aw+0.5*rho*CD*Ap*abs(Vw-V(n+1))*(Vw-V(n+1)))* (V(n+1));
Pa(n+1) = Power;
end
%---------------------------------------------------------------------
% the code below is to use the Euler method
% Xe and Ve are the displacement and the velocity from the Euler Method
Xe(1:ndt+1) = 0; % save displacement X from step 0 to step ndt
Ve(1:ndt+1) = 0;
PaE(1:ndt+1) = 0;
for n=1:ndt
% to calculate kx1 and kv1
ta = time(n);
aw = (omega *Um*cos(omega *ta))-((2*omega *Um/3)*cos(2*omega*ta)); % acceleration of the fluid
Vw = Um*sin(omega*ta)-(Um/3)*sin(2*omega*ta); % velocity of the fluid
t1 = CA*md*aw;
t2 = 0.5*rho*CD*Ap*abs(Vw-Ve(n))*(Vw-Ve(n));
t3 = -K*Xe(n)-c*Ve(n);
Xe(n+1) = Xe(n)+dt*Ve(n);
Ve(n+1) = Ve(n)+dt*(t1+t2+t3)/(m+CA*md);
PowerE =(CA*md*aw+0.5*rho*CD*Ap*abs(Vw-Ve(n+1))*(Vw-Ve(n+1)))* (Ve(n+1));
PaE(n+1) = PowerE;
end
%average of power
Pam(1:ndt+1) = 0;
PamE(1:ndt+1) = 0;
for n = 1:ndt
Pam(n) = mean (Pa);
PamE(n) = mean (PaE);
end
% write the results in to a file output.txt
fileID = fopen('output.txt','w');
for n=1:ndt+1
fprintf(fileID,'%12.6e %12.6e %12.6e %12.6e %12.6e\r\n', ...
time(n),X(n),V(n),Xe(n),Ve(n));
end
fclose(fileID);
% draw the displacement of the cylinder
figure (01);
plot(time(1:ndt),X(1:ndt),'-r');
xlabel('time (s)'); %honrizontal axis's label is 'x'
ylabel('X (m)'); %vertical axis's label is 'T'
hold;
plot(time(1:ndt),Xe(1:ndt),'-b');
legend ('RK method', 'Euler method');
% draw the velocity of the cylinder
figure (02)
plot(time(1:ndt),V(1:ndt),'-r');
xlabel('time (s)'); %honrizontal axis's label is 'x'
ylabel('V (m/s)'); %vertical axis's label is 'T'
hold;
plot(time(1:ndt),Ve(1:ndt),'-b');
legend ('RK method', 'Euler method');
figure (03)
plot(time(1:ndt),Pa(1:ndt),'-r');
xlabel('time (s)'); %honrizontal axis's label is 'x'
ylabel('Power (W)'); %vertical axis's label is 'T'
hold;
plot(time(1:ndt),PaE(1:ndt),'-b');
legend ('RK method', 'Euler method');
figure (04)
plot(time(1:ndt),Pam(1:ndt),'-r');
xlabel('time (s)'); %honrizontal axis's label is 'x'
ylabel('power(W)'); %vertical axis's label is 'T'
hold;
plot(time(1:ndt),PamE(1:ndt),'-b');
legend ('RK method', 'Euler method');
Vatsal Dhamelia
Vatsal Dhamelia 2020년 5월 11일
can you please help me fo plot KC vs Pam

댓글을 달려면 로그인하십시오.


Star Strider
Star Strider 2018년 3월 6일
One possibility:
x = linspace(2, 20, 25);
y = 8*x + 10;
figure(1)
plot(x, y, '-pg')
grid
xlabel('x')
ylabel('y')
title('Moody Chart')
See the documentation for the various functions for details on their uses.

카테고리

Help CenterFile Exchange에서 Fluid Dynamics에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by