Clear Filters
Clear Filters

Plotting a function given a range of inputs

53 views (last 30 days)
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.

Answers (2)

Elias Gule
Elias Gule on 6 Mar 2018
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 Comments
Vatsal Dhamelia
Vatsal Dhamelia on 11 May 2020
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 on 11 May 2020
can you please help me fo plot KC vs Pam

Sign in to comment.


Star Strider
Star Strider on 6 Mar 2018
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.

Categories

Find more on Fluid Dynamics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!