Plotting a function given a range of inputs

136 visualizzazioni (ultimi 30 giorni)
Jesse Green
Jesse Green il 6 Mar 2018
Commentato: Vatsal Dhamelia il 11 Mag 2020
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.

Risposte (2)

Elias Gule
Elias Gule il 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 Commenti
Vatsal Dhamelia
Vatsal Dhamelia il 11 Mag 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 il 11 Mag 2020
can you please help me fo plot KC vs Pam

Accedi per commentare.


Star Strider
Star Strider il 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.

Categorie

Scopri di più su Fluid Dynamics in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by