Plotting in 3D, based on beta and alpha.

Dear all,

I am trying to learn matlab and use it for my project.

In my code below I have made a plot P_o1 dependend from variable beta. I would like to make a 3D plot making it dependend from alpha aswell. Alpha is now set on 80 but I would like to plot it from 0 to 90 degrees. Can anyone help me and explain how to make such a plot? 🙏

Kind regards,

R

clear all
close all
clc


%basic assumption/variables
g       = 9.81;                 %Gravity constant 
d_50    = 0.1;                  %Rock/particle diameter(50% of the rocks) [m]
rho_w   = 1025;                 %Water density [kg/m^3]
rho_s   =2600;                  %Stone density [kg/m^3]      
rho_r   =(rho_s-rho_w)/rho_w;   %relative density [kg/m^3] also desribed as R_sd
vis_l   =0.0000013 ;            %Miedema model excel[m^2/s])
C_fr    =0.416;                 %assumption friction factor bed and pipe, miedema book p451 regime 3
w_a     =0;                     %added water to pipe[m^3/s]
P_i     =1000/3600;             %Amount of rock added to fallpipe m^3/sec  
alpha2  =w_a/P_i;               %Ratio water added to fallpipe  Lecture  notes C. van Rhee. 
Cd      =0.4;                   %Form coefficient of rock, Notes rhee p.12

%situation dependent variable/inputs
L       =20 ;                   %Length of pipe
D_p     = 1 ;                   %Pipe Diameter
lampda_w=d_50/D_p   ;           %Ratio particle size and pipe diameter--> used for hindered settling 
D_H     = D_p ;                 %Hydraulic diameter, assumde to be te same as Pipe diamter in this case --> cirular pipe
R_e     =(3*D_H)/vis_l;         %Reynolds number
alpha  = 80 ;                   %Angle pipe, 90 is vertical & 0 is Horizontal
C_vb    =0.4 ;                  %Volumetric bed concentration
beta    =pi/4;                  %Angle bed in pipe, Max input angle is pi/2!! described in figure 1, paper OMAW2014
Ap      =(pi/4.5)*D_p^2;        %Surface cross section pipe, described in figure 1, paper OMAW2014
A2      =Ap*(beta-sin(beta)*cos(beta)) ;    %Surface cross section of sliding bed, described in figure 1, paper OMAW2014
A1      =Ap-A2 ;                %Surface cross section of water in pipe, described in figure 1, paper OMAW2014
c       =A2/A1 ;                %Ratio surface cross section of water and sliding deck. 
v_l     =w_a/A1;                %Avarage velocity water in pipe[-], assumed that only A1 will be accerelated by adding the water. 
O_1     =D_p*(pi-beta);         %Length of water or the suspension in contact with the pipe eq 2 paper OMOW2014
O_2     =D_p*beta;              %Length of the fixed or sliding bed in contact with wall eq 3 paper OMOW2014
O_12    =D_p*sin(beta);         %Surface of bed in contact with water eq 4 paper OMOW2014
a_wilson=2.75;                  %Wilson factor p3 at eq 18 paper OMOW2014
n       =2.4 ;                  %The value chosen for the hindered sttling exponent n = 2.4, which is a value for high particle Reynolds numbers. This value is chosen because relative large particles are used for subsea rock installation (D = 0.02 - 0.1 [m]).


%These 'parts' are only used to prevent typos in the equation for P_o1 below. 

% part_1=rho_w*g*L*rho_r*((D_p^2)/4)*(beta-sin(beta)*cos(beta))*sind(alpha);
% part_2=(C_fr*rho_w*g*rho_r*C_vb*Ap*(beta-sin(beta)*cos(beta)))/(beta*D_p*pi);
% part_3=(a_wilson*1.325)/((log((0.27*d_50)/D_H+5.75/(R_e^0.9)))^2);
% part_4=(1/8)*rho_w*L*D_p*sin(beta);
% part_5=(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*(1-C_vb);
% part_6=D_p-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta));
% part_7=(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*C_vb;
% P_o=((sqrt((part_1-part_2)/(part_3*part_4)))/(1-part_5/part_6))*part_7 ;

figure (1)

betaaa =linspace(0,0.3*pi);
for i = 1:length(betaaa);
    beta = betaaa(i) ;
P_o1=((sqrt((rho_w*g*L*rho_r*((D_p^2)/4)*(beta-sin(beta)*cos(beta))*sind(alpha)-(C_fr*rho_w*g*rho_r*C_vb*Ap*(beta-sin(beta)*cos(beta)))/(beta*D_p*pi))/((a_wilson*1.325)/((log((0.27*d_50)/D_H+5.75/(R_e^0.9)))^2)*(1/8)*rho_w*L*D_p*sin(beta))))/(1-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*(1-C_vb)/(D_p-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta)))))*(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*C_vb ;
P_o1_list(i) = P_o1;
end
plot(betaaa,P_o1_list)
title('Production fallpipe assuming sliding bed, fixed angle')
xlabel('Beta size of sliding bed') 
ylabel('production')
grid on

Answer 1

[KSSV]

I strictly advice you to check your formul for P_o1.

%basic assumption/variables
g       = 9.81;                 %Gravity constant
d_50    = 0.1;                  %Rock/particle diameter(50% of the rocks) [m]
rho_w   = 1025;                 %Water density [kg/m^3]
rho_s   =2600;                  %Stone density [kg/m^3]
rho_r   =(rho_s-rho_w)/rho_w;   %relative density [kg/m^3] also desribed as R_sd
vis_l   =0.0000013 ;            %Miedema model excel[m^2/s])
C_fr    =0.416;                 %assumption friction factor bed and pipe, miedema book p451 regime 3
w_a     =0;                     %added water to pipe[m^3/s]
P_i     =1000/3600;             %Amount of rock added to fallpipe m^3/sec
alpha2  =w_a/P_i;               %Ratio water added to fallpipe  Lecture  notes C. van Rhee.
Cd      =0.4;                   %Form coefficient of rock, Notes rhee p.12

%situation dependent variable/inputs
L       =20 ;                   %Length of pipe
D_p     = 1 ;                   %Pipe Diameter
lampda_w=d_50/D_p   ;           %Ratio particle size and pipe diameter--> used for hindered settling
D_H     = D_p ;                 %Hydraulic diameter, assumde to be te same as Pipe diamter in this case --> cirular pipe
R_e     =(3*D_H)/vis_l;         %Reynolds number
alpha  = linspace(0,90)  ;                   %Angle pipe, 90 is vertical & 0 is Horizontal
C_vb    =0.4 ;                  %Volumetric bed concentration
beta    =pi/4;                  %Angle bed in pipe, Max input angle is pi/2!! described in figure 1, paper OMAW2014
Ap      =(pi/4.5)*D_p^2;        %Surface cross section pipe, described in figure 1, paper OMAW2014
A2      =Ap*(beta-sin(beta)*cos(beta)) ;    %Surface cross section of sliding bed, described in figure 1, paper OMAW2014
A1      =Ap-A2 ;                %Surface cross section of water in pipe, described in figure 1, paper OMAW2014
c       =A2/A1 ;                %Ratio surface cross section of water and sliding deck.
v_l     =w_a/A1;                %Avarage velocity water in pipe[-], assumed that only A1 will be accerelated by adding the water.
O_1     =D_p*(pi-beta);         %Length of water or the suspension in contact with the pipe eq 2 paper OMOW2014
O_2     =D_p*beta;              %Length of the fixed or sliding bed in contact with wall eq 3 paper OMOW2014
O_12    =D_p*sin(beta);         %Surface of bed in contact with water eq 4 paper OMOW2014
a_wilson=2.75;                  %Wilson factor p3 at eq 18 paper OMOW2014
n       =2.4 ;                  %The value chosen for the hindered sttling exponent n = 2.4, which is a value for high particle Reynolds numbers. This value is chosen because relative large particles are used for subsea rock installation (D = 0.02 - 0.1 [m]).

betaaa =linspace(0,0.3*pi);
[alpha,beta] = meshgrid(alpha,betaaa) ; 

figure (1)
P_o1=((sqrt((rho_w*g*L*rho_r*((D_p^2)/4).*(beta-sin(beta).*cos(beta)).*sind(alpha).....
      -(C_fr*rho_w*g*rho_r*C_vb*Ap*(beta-sin(beta)*cos(beta)))./(beta*D_p*pi))./((a_wilson*1.325)./((log((0.27*d_50)./D_H+5.75/(R_e^0.9)))^2)*(1/8)*rho_w*L*D_p*sin(beta))))./(1-(pi/4)*D_p^2*(beta-sin(beta).*cos(beta))*(1-C_vb)./(D_p-(pi/4)*D_p^2*(beta-sin(beta).*cos(beta))))).*(pi/4)*D_p^2.*(beta-sin(beta)*cos(beta))*C_vb ;
    

surf(alpha,beta,P_o1)
title('Production fallpipe assuming sliding bed, fixed angle')
xlabel('Beta size of sliding bed')
ylabel('production')
grid on

[R]

Thank you for your help, the graph you plotted is excatly something I am looking for. However in none of my variables(beta, alpha and p_o1) for the P_o1 equation can be negative, so there might be something wrong with your plot, I could not find what it was.

I did find a minor error in my P_o1 equation, however my output did not change that much. I deleted all the unneccesairy lines in my code and made the following:

Now I plotted my P_o1 six times with different alpha values, all the values seem to be in the right order of magnitude but I would like to combine them in one 3D plot, do you know how to do that? Now I have alpha 10/20/30/40/50/60, eventually I need to have them for all angles from +/- 10 to +/- 70.

clear all
close all
clc


%basic assumption/variables
g       = 9.81;                 %Gravity constant 
d_50    = 0.1;                  %Rock/particle diameter(50% of the rocks) [m]
rho_w   = 1025;                 %Water density [kg/m^3]
rho_s   =2600;                  %Stone density [kg/m^3]      
rho_r   =(rho_s-rho_w)/rho_w;   %relative density [kg/m^3] also desribed as R_sd
vis_l   =0.0000013 ;            %Miedema model excel[m^2/s])
C_fr    =0.416;                 %assumption friction factor bed and pipe, miedema book p451 regime 3
w_a     =0;                     %added water to pipe[m^3/s]
P_i     =1000/3600;             %Amount of rock added to fallpipe m^3/sec  
Cd      =0.4;                   %Form coefficient of rock, Notes rhee p.12

%situation dependent variable/inputs
L       =20 ;                   %Length of pipe
D_p     = 1.2 ;                 %Pipe Diameter
lampda_w=d_50/D_p   ;           %Ratio particle size and pipe diameter--> used for hindered settling 
D_H     = D_p ;                 %Hydraulic diameter, assumde to be te same as Pipe diamter in this case --> cirular pipe
R_e     =(3*D_H)/vis_l;         %Reynolds number
C_vb    = 0.4 ;                 %Volumetric bed concentration
beta    =pi/4;                  %Angle bed in pipe, Max input angle is pi/2!! described in figure 1, paper OMAW2014
Ap      =(pi/4)*D_p^2;          %Surface cross section pipe, described in figure 1, paper OMAW2014
a_wilson=2.75;                  %Wilson factor p3 at eq 18 paper OMOW2014
n       =2.4 ;                  %The value chosen for the hindered sttling exponent n = 2.4, which is a value for high particle Reynolds numbers. This value is chosen because relative large particles are used for subsea rock installation (D = 0.02 - 0.1 [m]).



alpha1  = 10 ;                   %Angle pipe, 90 is vertical & 0 is Horizontal
figure(1)
betaaa =linspace(0.01,0.2*pi);
for i = 1:length(betaaa);
    beta = betaaa(i) ;
P_o1=((sqrt((rho_w*g*L*rho_r*((D_p^2)/4)*(beta-sin(beta)*cos(beta))*sind(alpha1)-(C_fr*rho_w*g*cosd(alpha1)*rho_r*C_vb*Ap*(beta-sin(beta)*cos(beta)))/(beta*D_p*pi))/((a_wilson*1.325)/((log((0.27*d_50)/D_H+5.75/(R_e^0.9)))^2)*(1/8)*rho_w*L*D_p*sin(beta))))/(1-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*(1-C_vb)/(D_p-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta)))))*(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*C_vb ;
P_o1_list(i) = P_o1;
end
plot(betaaa,P_o1_list)
title('Production fallpipe assuming sliding bed, fixed angle')
xlabel('Beta size of sliding bed') 
ylabel('production')
grid on


alpha2  = 20 ;                   %Angle pipe, 90 is vertical & 0 is Horizontal
figure(2)
betaaa =linspace(0.01,0.2*pi);
for i = 1:length(betaaa);
    beta = betaaa(i) ;
P_o1=((sqrt((rho_w*g*L*rho_r*((D_p^2)/4)*(beta-sin(beta)*cos(beta))*sind(alpha2)-(C_fr*rho_w*g*cosd(alpha2)*rho_r*C_vb*Ap*(beta-sin(beta)*cos(beta)))/(beta*D_p*pi))/((a_wilson*1.325)/((log((0.27*d_50)/D_H+5.75/(R_e^0.9)))^2)*(1/8)*rho_w*L*D_p*sin(beta))))/(1-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*(1-C_vb)/(D_p-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta)))))*(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*C_vb ;
P_o1_list(i) = P_o1;
end
plot(betaaa,P_o1_list)
title('Production fallpipe assuming sliding bed, fixed angle')
xlabel('Beta size of sliding bed') 
ylabel('production')
grid on

alpha3  = 30 ;                   %Angle pipe, 90 is vertical & 0 is Horizontal
figure(3)
betaaa =linspace(0.01,0.2*pi);
for i = 1:length(betaaa);
    beta = betaaa(i) ;
P_o1=((sqrt((rho_w*g*L*rho_r*((D_p^2)/4)*(beta-sin(beta)*cos(beta))*sind(alpha3)-(C_fr*rho_w*g*cosd(alpha3)*rho_r*C_vb*Ap*(beta-sin(beta)*cos(beta)))/(beta*D_p*pi))/((a_wilson*1.325)/((log((0.27*d_50)/D_H+5.75/(R_e^0.9)))^2)*(1/8)*rho_w*L*D_p*sin(beta))))/(1-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*(1-C_vb)/(D_p-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta)))))*(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*C_vb ;
P_o1_list(i) = P_o1;
end
plot(betaaa,P_o1_list)
title('Production fallpipe assuming sliding bed, fixed angle')
xlabel('Beta size of sliding bed') 
ylabel('production')
grid on


alpha4  = 40 ;                   %Angle pipe, 90 is vertical & 0 is Horizontal
figure(4)
betaaa =linspace(0.01,0.2*pi);
for i = 1:length(betaaa);
    beta = betaaa(i) ;
P_o1=((sqrt((rho_w*g*L*rho_r*((D_p^2)/4)*(beta-sin(beta)*cos(beta))*sind(alpha4)-(C_fr*rho_w*g*cosd(alpha4)*rho_r*C_vb*Ap*(beta-sin(beta)*cos(beta)))/(beta*D_p*pi))/((a_wilson*1.325)/((log((0.27*d_50)/D_H+5.75/(R_e^0.9)))^2)*(1/8)*rho_w*L*D_p*sin(beta))))/(1-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*(1-C_vb)/(D_p-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta)))))*(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*C_vb ;
P_o1_list(i) = P_o1;
end
plot(betaaa,P_o1_list)
title('Production fallpipe assuming sliding bed, fixed angle')
xlabel('Beta size of sliding bed') 
ylabel('production')
grid on

alpha5  = 50 ;                   %Angle pipe, 90 is vertical & 0 is Horizontal
figure(5)
betaaa =linspace(0.01,0.2*pi);
for i = 1:length(betaaa);
    beta = betaaa(i) ;
P_o1=((sqrt((rho_w*g*L*rho_r*((D_p^2)/4)*(beta-sin(beta)*cos(beta))*sind(alpha4)-(C_fr*rho_w*g*cosd(alpha4)*rho_r*C_vb*Ap*(beta-sin(beta)*cos(beta)))/(beta*D_p*pi))/((a_wilson*1.325)/((log((0.27*d_50)/D_H+5.75/(R_e^0.9)))^2)*(1/8)*rho_w*L*D_p*sin(beta))))/(1-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*(1-C_vb)/(D_p-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta)))))*(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*C_vb ;
P_o1_list(i) = P_o1;
end
plot(betaaa,P_o1_list)
title('Production fallpipe assuming sliding bed, fixed angle')
xlabel('Beta size of sliding bed') 
ylabel('production')
grid on

alpha6  = 60 ;                   %Angle pipe, 90 is vertical & 0 is Horizontal
figure(6)
betaaa =linspace(0.01,0.2*pi);
for i = 1:length(betaaa);
    beta = betaaa(i) ;
P_o1=((sqrt((rho_w*g*L*rho_r*((D_p^2)/4)*(beta-sin(beta)*cos(beta))*sind(alpha4)-(C_fr*rho_w*g*cosd(alpha4)*rho_r*C_vb*Ap*(beta-sin(beta)*cos(beta)))/(beta*D_p*pi))/((a_wilson*1.325)/((log((0.27*d_50)/D_H+5.75/(R_e^0.9)))^2)*(1/8)*rho_w*L*D_p*sin(beta))))/(1-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*(1-C_vb)/(D_p-(pi/4)*D_p^2*(beta-sin(beta)*cos(beta)))))*(pi/4)*D_p^2*(beta-sin(beta)*cos(beta))*C_vb ;
P_o1_list(i) = P_o1;
end
plot(betaaa,P_o1_list)
title('Production fallpipe assuming sliding bed, fixed angle')
xlabel('Beta size of sliding bed') 
ylabel('production')
grid on